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Abstract 


Solutions to ever-larger structural optimization problems are desired. 
However, computational resources are strained to meet this need. New meth- 
ods will be required to solve ever larger problems. The present approaches 
to solving large-scale problems involve approximations for the constraints of 
structural optimization problems and/or decomposition of the problem into 
multiple subproblems that can be solved in parallel. An area of game theory, 
equilibrium programming (also know as non-cooperative game theory), can be 
used to unify these existing approaches from a theoretical point of view (consid- 
ering the existence and optimality of solutions), and be used as a framework for 
the development of new methods for solving large-scale optimization problems. 
Equilibrium programming theory is described, and existing design techniques 
such as fully stressed design and constraint approximations are shown to fit 
within its framework. Two new structural design formulations are also derived. 
The first new formulation is another approximation technique wdiich is a gen- 
eral updating scheme for the sensitivity derivatives of design constraints. The 
second new formulation uses a substructure-based decomposition of the struc- 
ture for analysis and sensitivity calculations. The new formulations are utilized 
for problems ranging from a simple ten-bar truss to a 348 design- variable opti- 
mization of a high-speed civil transport. Significant computational benefits of 
the new formulations compared with a conventional method are demonstrated. 


1 



Acknowledgements 


This paper is essentially a dissertation submitted to the faculty of the School 
of Engineering and Applied Science of The George Washington University in 
partial satisfaction of the requirements for the degree of Doctor of Science. I 
wish to express my sincere gratitude to my advisors, Dr. Raphael Haftka, and 
Dr. John Whitesides. I also would like to thank my family, especially my 
wife Joan, and children Brian, Sarah, Shea, and Samantha for their patience 
during my long tenure as a graduate student, and for their acceptance of the 
demands it placed on my time. Mr. Carl Martin was of great assistance 
in development of the SLP runstream used in the example problems, and in 
setting up and running some of the example problems. The contribution of 
the Boeing Company in providing the vehicle example problem is appreciated. 
The support and encouragement of the NASA Langley Research Center, in 
particular my branch head Dr. Charles Camarda and my former branch heads 
Dr. James Starnes and Dr. Donald Rummler, is gratefully acknowledged. 



Table of Contents 


Abstract i 

Acknowledgements ii 

Table of Contents iii 

List of Tables v 

List of Figures vi 

Nomenclature viii 

Chapter 1. Introduction 1 

1.1 Background on Structural Synthesis 1 

1.2 Scope of Present Study 8 

Chapter 2. Equilibrium Programming Background 13 

2.1 Equilibrium Programming Background 13 

2.2 Features of Equilibrium Programming Problems 15 

2.3 Existence of an Equilibrium Point 16 

2.4 Some Equilibrium Point Properties 18 

2.5 Solution Methodologies 20 

Chapter 3. Development of Equilibrium Programming Structural- 
Design Formulations 22 

3.1 Optimal Structural-Sizing Problem Statement 24 

3.2 Structural-Design Formulation 1 26 

3.3 Structural- Design Formulation 2 31 

3.4 Structural- Design Formulation 3 34 

3.4.1 Structural- Design Formulation Using Updated Sensitivity Deriva- 
tives 36 

iii 


3.5 Design Formulation 4 38 

3.5.1 Derivation of the Substructure Analysis Method 40 

3.5.2 Derivation of EP Design Formulation 4 47 

Chapter 4. Numerical Results for Equilibrium Programming Structural- 
Design Formulations 57 

4.1 Implementation of Numerical Algorithms 57 

4.1.1 General Implementation Issues 58 

4.1.2 Implementation Issues for Design Formulation 3 61 

4.1.3 Implementation Issues for Design Formulation 4 63 

4.2 Description of the Example Problems 64 

4.2.1 Ten-Bar-Truss Example Problem 65 

4.2.2 High-Speed Civil Transport Wing Example Problem 66 

4.2.3 High-Speed Transport Vehicle Example Problem 68 

4.2.4 Transmission Tower Example Problem 69 

4.3 Results for the Example Problems 71 

4.3.1 Ten-Bar- Truss Example Problem Results 71 

4.3.2 High-Speed Civil Transport Wing Example Problem Results . 77 

4.3.3 High-Speed Transport Vehicle Example Problem Results ... 82 

4.3.4 Transmission Tower Example Problem Results 87 

Chapter 5. Concluding Remarks 91 

Bibliography 95 


iv 



List of Tables 


4.1 Comparison of normalized CPU timing results for the trans- 
mission tower example using EP structural-design formulation 
4 with two and four substructures, and using a conventional 
SLP approach for the entire structure 89 


v 



List of Figures 


3.1 Design space for a rod constrained between two rigid walls un- 
dergoing a temperature change 31 

3.2 Example wing structural finite element model decomposed into 

two substructures 42 

4.1 Flowchart of the solution procedure for structural design formu- 
lation 3 62 

4.2 Schematic of geometry, loads, and displacement limits for the 

ten-bar-truss example problem 66 

4.3 Finite element model of high-speed civil transport wing with 

upper cover panels removed 67 

4.4 Half-symmetric, finite element model of the high-speed transport 

vehicle 69 

4.5 Schematic of the transmission tower example geometry, loading, 

and division into substructures 70 

4.6 Comparison of penalized objective function iteration histories for 

the ten-bar truss with different initial move-limit-control factors 
using EP structural-design formulation 3, and using a conven- 
tional SLP approach 73 

4.7 Comparison of objective function iteration histories for the ten- 

bar truss with different initial move-limit-control factors using 
EP structural-design formulation 3, and using a conventional 
SLP approach 75 

4.8 Comparison of iteration histories for displacement constraint on 

<5i for the ten-bar truss with different initial move-limit-control 
factors using EP structural-design formulation 3, and using a 
conventional SLP approach 75 

4.9 Comparison of penalized objective function iteration histories 

using a CPU time ordinate for the ten-bar truss with different 
initial move-limit-control factors using EP structural-design for- 
mulation 3, and using a conventional SLP approach 76 

4.10 Comparison of penalized objective function iteration history for 
the high-speed civil transport wing with different initial move- 
limit-control factors using EP structural-design formulation 3, 

and using a conventional SLP approach 78 


vi 



4.11 Objective function and wing tip displacement constraint itera- 

tion histories for the high-speed civil transport wing using EP 
structural-design formulation 3, and using a conventional SLP 
approach 81 

4.12 Comparison of penalized objective function iteration history us- 

ing a CPU time ordinate for the high-speed civil transport wing 
with different initial move-limit-control factors using EP structural- 
design formulation 3, and using a conventional SLP approach. . 81 

4.13 Comparison of penalized objective function iteration history for 
the high-speed transport vehicle using EP structural-design for- 
mulation 3, using modified EP structural-design formulation 3, 

and using a conventional SLP approach 83 

4.14 Comparison of objective function and most critical stress con- 

straint iteration history for the high-speed transport vehicle us- 
ing EP structural-design formulation 3, using modified EP structural- 
design formulation 3, and using a conventional SLP approach. . 86 

4.15 Comparison of penalized objective function iteration history us- 

ing a CPU time ordinate for the high-speed transport vehicle 
using EP structural-design formulation 3, using modified EP 
structural-design formulation 3, and using a conventional SLP 
approach 86 

4.16 Comparison of objective function and most critical stress con- 

straint iteration history for the transmission tower example us- 
ing EP structural-design formulation 4 with two and four sub- 
structures, and using a conventional SLP approach for the entire 
structure 88 


vii 



Nomenclature 


Column vectors are generally denoted by lower case symbols typed in boldface. 
The notation df/da for an arbitrary scalar / and a vector a = (ai, . . . , a n ) T 
denotes the gradient of scalar / with respect to vector a. This gradient is 
expressed as a row vector (i.e., df/da = ( df jda \ , . . . , df/da n )). The notation 
di/da for arbitrary vectors a and f = (/i,-..,/ m ) T denotes the gradient of 
vector f with respect to vector a. This gradient is expressed as an m x n 
matrix in which the component in column j of row i is given by dfi/daj. The 
notation fij denotes the j^ 1 component of the vector f*. 

a update vector 

A cross-sectional area, sometimes (in formulation 4) the set 
of 4-tuples used to define compatibility of substructures 

A vector of cross-sectional areas 

B j matrices used to define compatible displacements at sub- 
structure interfaces 

g vector of inequality constraint functions 

ej unit vector with length of vector v and having unity for 
the j’th entry 

/ objective function in nonlinear programming 

F vector of external nodal forces 

h vector of equality constraint functions 

k constant of proportionality 

K penalty coefficient used in penalized objective function P 

viii 



K, K global stiffness matrix for a structural-response prob- 
lem, augmented global stiffness matrix for a structural- 
response subproblem (see equations (3.28) and (3.34)) 

m number of substructures that would have rigid-body mo- 

tions if separated from the rest of the structure, used in 
formulation 4 

M number of EP subproblems, sometimes number of sub- 

structures in formulation 4 

M matrix defining the substructure coordination problem 

(see equations (3.29) and (3.36)) 

n number of displacement degrees of freedom 

N vector of stress resultants within every element of the 

structure 

P penalized objective function (see equation (4.2)) 

R matrix containing substructure rigid-body modes 

u vector of nodal displacements from a structural- response 

problem 

U vector of nodal displacements orthogonal to rigid-body 

modes 

v vector of design variables for nonlinear programming, or 

vector of structural-sizing design variables in equilibrium 
programming structural-sizing subproblems 

Wk weights used in update subproblem (see statement (3.12)) 

W weight of a structure 

x, Xj, Xj vectors of all design variables for equilibrium program- 
ming, design variables of subproblem i , and all design 
variables except those from subproblem i , respectively 

a vector of rigid-body mode amplitudes 

6 displacement value 


A/, Avk difference quantities used in update equation (3.13) 
AT temperature rise 
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A vector of Lagrange multipliers for inequality constraints, 

sometimes Lagrange multipliers for constraint on com- 
patible substructure interface displacements 

jj, vector of Lagrange multipliers for equality constraints, 

sometimes vector of Lagrange multipliers for constraint 
on orthogonality of u nodal displacements with rigid- 
body modes 

<r vector of stresses within every element of a structure 

Subscripts and Subscripts: 

A denotes an approximation of a function 

i denotes i’th equilibrium programming subproblem 

I denotes interface quantities in design formulation 4 

l, u denotes lower and upper bounds, used to denote move 

limits 

max maximum 

min minimum 

new, old new and old values, respectively; used in defining update 
subproblem (3.12) 

N denotes functional dependence on stress resultants 

s side constraints 

T transpose of a vector or matrix 

u denotes displacement constraints, or functional depen- 

dence on displacements 

<r denotes stress and local buckling constraints, or func- 

tional dependence on stresses 

* denotes an equilibrium programming solution (i.e., the 

equilibrium point), an optimal solution, active con- 
straints, or Lagrange multipliers corresponding to active 
constraints 
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Chapter 1 


Introduction 


Solutions to ever-larger structural optimization problems are needed as struc- 
tural optimization procedures are applied to solve practical design problems 
and the analysis models for these structures grow increasingly detailed and 
complex. Current applications of structural optimization to large-scale design 
problems require significant computational resources. New structural optimiza- 
tion methods are required to solve ever-larger problems. The use of a gener- 
alized theoretical framework for the development of structural optimization 
methods is investigated in the present dissertation. 

1.1 Background on Structural Synthesis 

Nonlinear mathematical programming (NLP) is now extensively used for a type 
of optimal structural design that is called structural sizing. The solution to a 
NLP-based, structural-sizing design problem minimizes an objective function, 
such as weight, and satisfies a set of design constraints, such as minimum gauge 
and local stress limits. This solution is obtained varying a set of structural- 
sizing design variables, such as skin thicknesses and stiffener heights, within 
their allowable ranges. The original usage of NLP for structural sizing coupled 
with finite element structural analyses was due to Schmit in 1960 (see [50] for a 
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description of the historical developments). Schmit called his resulting method 
structural synthesis. The structural synthesis method had several advantages 
over the state of the art at that time (i.e., fully stressed design methods and 
simultaneous failure mode methods). These advantages include such features 
as: 1) limits on member sizes were easily incorporated; 2) a true optimum design 
could be found subject to multiple load cases; 3) numerous failure modes could 
be considered simultaneously; and 4) objective functions other than weight 
were possible. The disadvantage of the method was the large demand placed 
on computational resources. 

In the original structural synthesis approach, the structural analysis re- 
quired to compute the design constraints was subordinate to the optimization 
algorithm, and a structural analysis was performed for every change in the 
values of the structural-sizing design variables in the design process. Thus, 
the process was computationally expensive. Several methods for overcoming 
this disadvantage such as the use of approximation techniques, reduction tech- 
niques, and decomposition techniques were presaged as early as 1971 by Pope 
and Schmit [40]. 

Approximation techniques are utilized to loosen the coupling of the 
structural analysis from the structural optimization procedure. Computation- 
ally efficient, approximate models for the dependence of the structural response 
and the constraint functions on the design variables are developed from the 
structural analysis and the sensitivity derivatives of the response. These ap- 
proximate models are used in place of the full structural analysis and constraint 
evaluation in the structural optimization procedure. After the structural op- 
timization procedure converges, the approximate models for the structural re- 
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sponse are reformulated using the resulting set of design variables, and the cycle 
is repeated until the variation of results between cycles is below a predefined 
tolerance. The accuracy of using linear Taylor series approximate analyses was 
evaluated by Storaasli and Sobieszczanski in [63] using a finite element model 
of a subsonic transport fuselage section. However, Schmit and Farshi defined 
the approximation concepts commonly used today in [45], and the further de- 
velopment of these concepts by Schmit and Miura [48] significantly improved 
the computational efficiency of NLP structural-sizing design methods. Many of 
the numerous approximation techniques that have been developed are summa- 
rized in [2]. Thus, approximation-based NLP structural-sizing design methods 
are now well developed, and are implemented in several commercial structural 
analysis codes. 

Reduction techniques are utilized to reduce the size of a structural op- 
timization problem to make it more tractable. In one approach, the original 
set of design variables is replaced by a sum of fixed basis vectors for the design 
variables (chosen by the user) multiplied by undetermined coefficients. If the 
number of basis vectors is smaller than their length, these coefficients form a 
new set of design variables that is smaller than the original set of design vari- 
ables (see [39]). However, the best choice of basis vectors is not known a priori , 
and the quality of the resulting design depends on the user’s skill in choosing 
good basis vectors. In another reduction technique, the constraint functions 
are combined into one or more cumulative constraint functions. This approach 
can significantly reduce the number of constraints. Cumulative constraints may 
make the adjoint sensitivity method more attractive by reducing the number of 
constraints to a number significantly below the number of design variables ([19], 



p. 267), and they are also utilized in some of the decomposition methods de- 
scribed subsequently. The Kreisselmeier-Steinhauser (K-S) function described 
in [27] is one example of a cumulative constraint function. 

Another approach that has been investigated to improve the computa- 
tional efficiency of structural sizing is the decomposition of the problem into a 
number of smaller subproblems that can be analyzed in parallel. An excellent 
summary of the different ways decomposition techniques are utilized to solve 
engineering design problems is given by Barthelemy [5] where decomposition 
methods are classified based on how the design variables are coupled through 
the problem constraints. This classification technique is similar to the method 
for classifying large linear programs ( [29], pp. 117 - 122). Decomposition tech- 
niques become computationally advantageous if the overall problem can be 
reformulated into subproblems having sets of design variables and constraints 
that are nearly disjoint. In practice, there is usually some coupling between 
the subproblems. If some design variables are common to the constraints of the 
different subproblems, they are called coupling variables. If some constraints 
depend on design variables from more than one subproblem, they are called 
coupling constraints. In linear programming, problems with coupling variables, 
but without coupling constraints, are called angular, and problems with cou- 
pling constraints, but no coupling variables, are called dual-angular. 

Some of the earliest descriptions of general processes for decomposition 
of large systems into subsystems are given by Mesarovick, et al. in [32] for cer- 
tain types of hierarchical systems. A hierarchical system can be decomposed 
into a tree or “family” of subsystems in which the “children” subsystems at a 
given level are independent of their “siblings” at the same level, and the cou- 
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pling between these sibling subsystems is through the “parent” subsystem. The 
parent subsystem controls its children subsystems using a set of coordination 
variables. Several methods of forming coordinable subproblems are described 
in [32]. These methods ensure that when a solution to the decomposed prob- 
lem is found, it is also a solution to the original problem. The hierarchical 
decomposition methods of Mesarovick become feasible when the design vari- 
ables of the problem can be grouped into sets that are only coupled through a 
few constraint equations (i.e., when there is constraint coupling). One common 
method for coordination when there is constraint coupling is that of goal coor- 
dination. In goal coordination, the constraints in a subproblem that have no 
explicit dependence on the sibling subsystems are treated explicitly, while the 
coupling constraint equations are multiplied by their Lagrange multipliers and 
are added to the objective functions, forming a partial Lagrangian function. 
The parent subsystem, which is the dual of the original problem, determines 
the optimum value of these Lagrange multipliers. Because goal coordination is 
essentially a primal-dual optimization method, the objective functions of the 
child subsystems must be convex, or must be “convexified” [7] for the method 
to yield satisfactory solutions [17]. This approach has been utilized for plastic 
design problems by Kaneko and Ha in [23]. 

The most common decomposition approach for structures treats a struc- 
ture in a hierarchical way with increasing levels of detail at each lower level 
in the hierarchy. For example, a discrete region - commonly referred to as a 
substructure - of a stiffened structure may be represented at an upper level in 
the hierarchy by smeared stiffnesses, while the region is modeled with discrete 
geometry at a lower level. The decomposition of a structure into substructures 
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to define optimization regions is described in this section. The decomposi- 
tion of a structure into substructures for analysis and optimization purposes is 
described in a subsequent chapter. 

Two early approaches that utilized two levels of representation for sub- 
structures and that performed searches in the lower-level design space are de- 
scribed by Giles in [16] for wing substructures, and by Kirsch, et al. in [25] for 
more general substructures. Schmit and Ramanathan [49] describe a two-level 
approach in which optimization is used within both levels. The upper-level 
subproblem utilizes weight as the objective function, and area (for trusses) or 
smeared orthotropic layer thicknesses (for stiffened skins) as design variables. 
The upper-level constraints are the side constraints (lower and upper limits) 
on the design variables, displacement constraints, system buckling constraints, 
and stress constraints. In the upper-level stress constraints, the compressive 
allowable stress is the maximum of a fixed lower stress limit and of the local 
buckling and crippling allowables that are determined from the lower-level de- 
sign of the previous iteration. The lower-level subproblems utilize the discrete 
geometry of substructures (such as stiffener blade height) as design variables, 
and minimize the difference of the structural stiffness determined using the 
lower-level design variables from the stiffness determined using the upper-level 
design variables. The constraints at this level are discrete geometry side con- 
straints, stress constraints, and local buckling and crippling constraints. At 
convergence of the two-level approach, the stiffness of the structure is the same 
whether it is determined by the upper-level or the lower-level design variables. 
A generalization of this approach that allows for laminated composite struc- 
tures is described by Schmit and Merhinfar in [47]. In these three two-level 
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approaches there is no coordination or direct coupling between the lower-level 
designs of substructures, and the resulting designs are not necessarily optimal. 

The linear decomposition method of Sobieski [58] is a widely investi- 
gated hierarchical scheme applied to structural optimization that may approach 
a locally optimal solution. (Although Kirsch [26] and Kirsch and Moses [24] 
developed rigorous, substructure-based approaches for hierarchical decomposi- 
tion of a structural optimization problem prior to the work of Sobieski, these 
methods are not finite-element based, and have only been applied to simple 
problems.) In the linear decomposition method, the cumulative constraint vi- 
olation for the lower-level constraints of Schmit and Ramanathan is minimized 
in each lower-level subproblem. The only explicit constraints in a lower-level 
subproblem are equality constraints that ensure that the stiffness determined 
using the lower-level design variables equals the structural stiffness determined 
using the upper-level design variables. An optimum sensitivity analysis with 
respect to the upper-level design variables, as described in references [4] and 
[55], is performed for these lower- level, cumulative-constraint, objective func- 
tions. At the upper level, the problem formulation is as described in Schmit 
and Ramanathan except that the upper-level stress constraints are replaced 
by linear approximations to the optimal cumulative constraints using the op- 
timization results and the optimal sensitivity analyses of the lower level. This 
method is demonstrated in [56], generalized for multiple levels in [57], and 
applied to a large transport aircraft sizing problem in [67]. The method is 
combined with approximation methods for the constraints in [3], and two more 
recent variations of the method are described in [61]. 

Since the linear decomposition method is essentially the decomposition 
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of a problem using coupling variables (the upper-level design variables enter 
into the lower- level subproblems as parameters), other solution methods that 
allow coupling variables can be used to form alternate solution strategies. The 
penalty approach of Haftka [21], and the method of Thareja and Haftka [64] 
are two such examples, which interestingly enough can be efficiently solved at a 
single level. More recently, nonhierarchical decomposition methods have been 
developed to simplify the multidisciplinary optimization process. These meth- 
ods utilize the notion of a global sensitivity analysis of interacting subsystems 
[60] to determine the total sensitivity of interdependent subsystems to changes 
in problem parameters, and concurrent subspace methods along with an over- 
all coordinating problem for optimization (see references [8], [54], and [59]). 
Global sensitivity analysis has also been shown to be useful for hierarchical 
structural analysis by Padula and Polignone in [38]. 

1.2 Scope of Present Study 

In the present dissertation, efficient methods for solving structural optimiza- 
tion problems are developed by approaching an optimization problem using a 
theoretical framework that is more general than the nonlinear programming 
theory that forms the basis of the structural optimization methods described 
in the previous section. The approach of using a different, and perhaps un- 
conventional, theoretical framework to take a fresh look at an old problem 
is common in mathematics. Utilizing a generalized theoretical framework in 
studying structural optimization allows for a new perspective of the existing 
solution methods and suggests new approaches that are more efficient than the 
existing approaches. 
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Equilibrium programming (EP), or non-cooperative game theory, is the 
generalized theoretical framework utilized in the present dissertation. A more 
detailed overview of EP is presented in the next chapter, but in brief, equilib- 
rium programming is a theory that describes the behavior of multiple, interact- 
ing systems that can each be described as NLP problems. An important advan- 
tage of utilizing an equilibrium programming framework for developing struc- 
tural optimization methods over ad-hoc approaches is that a well-developed 
theory concerning existence and optimality of solutions is available in the lit- 
erature. Thus, the development of the equilibrium programming structural 
design formulations in the present dissertation was guided by, and benefitted 
from, the theoretical structure provided by EP. 

Additional guidance for the development of computationally efficient, 
equilibrium programming structural design formulations is obtained by study- 
ing the reasons for the success of the commonly used design methods. As 
indicated in the previous section, the commonly used methods for solving large- 
scale structural optimization problems include approximations for the design 
constraints, various methods for reducing the size of the problem, and decom- 
position of the problem into multiple subproblems that can be solved in paral- 
lel. These methods have been very successful in improving the computational 
efficiency of finite-element-based structural analysis and optimal structural de- 
sign. A study of the computational advantages provided by these methods can 
suggest steps to further improve computational efficiency. 

Using approximation concepts, such as simple, approximate equations 
to describe the design constraint functions, the number of expensive finite el- 
ement analyses required in the NLP solution is reduced significantly. But a 
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primary reason for the improved efficiency is that the number of sensitivity 
analyses, which supply sensitivity derivatives for the gradient-based search 
methods commonly used, is dramatically decreased. It is often stated that 
the sensitivity derivatives obtained by the semi-analytic direct sensitivity anal- 
ysis method (initially developed in [14]) are inexpensive because the factors 
for the stiffness matrix in the displacement sensitivity equation are available 
from the structural analysis. However, the right-hand side of the displacement 
sensitivity equation has as many columns as there are design variables. As a 
result, the computation of the right-hand side for the semi-analytic sensitivity 
equation, the back substitution to determine the matrix of displacement sen- 
sitivity derivatives, and the use of this matrix in chain rule or finite difference 
calculations to determine stress sensitivity derivatives can be very expensive. 
In addition, the computational effort for computing these sensitivity derivatives 
increases in almost direct proportion to the number of load cases (however, as 
the number of load cases increases, the adjoint sensitivity method may become 
more computationally efficient than the direct method). The dominance of sen- 
sitivity analysis cost over the structural-response analysis cost is alluded to in 
[3], and it is explicitly demonstrated for the larger example problems described 
subsequently. Thus, as problems become more complex, a fundamental need is 
to reduce the number of sensitivity derivative calculations. 

An understanding of the primary advantages of decomposition of an 
optimization problem can also guide development of methods having increased 
computational efficiencies. A fundamental advantage of the linear decompo- 
sition method is that when a number of lower-level subproblems are to be 
solved, they can be solved in parallel. Another advantage, which is not often 
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highlighted, is that the displacement sensitivity derivatives calculated axe sen- 
sitivity derivatives with respect to a small set of upper-level design variables. 
In other words, these derivatives are with respect to a relatively small number 
of substructure stiffnesses, not with respect to the more numerous detailed ge- 
ometry parameters. Thus, the number of displacement sensitivity derivatives 
required can be greatly reduced using the linear decomposition method. Even 
with this advantage, it was determined for the example problem in [67] that 
36% of the total computer time for problem solution was utilized in the finite- 
element-based analysis and sensitivity derivative calculations that involved only 
five of the 1300 design variables! Thus, although the new optimization methods 
developed in the present dissertation are derived as equilibrium programming 
formulations, one of the fundamental approaches for the development of these 
methods is to find ways to reduce the number of sensitivity analyses required, 
and to reduce the size of the sensitivity analysis problems. 

In subsequent chapters, equilibrium programming theory is summarized, 
four equilibrium programming structural design formulations are developed, 
and numerical results for two of these formulations are presented. Linear struc- 
tural analysis is assumed in the derivations and test problems, and the con- 
straints considered are side constraints, displacement constraints, stress and 
local buckling constraints. Specifically, some of the fundamental properties of 
equilibrium programming are summarized in chapter 2. In chapter 3, the basic 
equations governing finite-element-based structural analysis and optimization 
are described, and two commonly used design methods, fully stressed design 
and constraint approximations, are developed as EP design formulations. Two 
new EP design formulations are also derived in the chapter. The first new EP 
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design formulation utilizes approximate, updated sensitivity derivatives that 
replace the constraint sensitivity derivatives determined by a traditional sen- 
sitivity analysis at a small fraction of the cost. The second new EP design 
formulation is a substructure-based decomposition method in which the sensi- 
tivity analysis subproblems are greatly reduced in size, and can be solved in 
parallel. The method of derivation of this formulation ensures the existence 
and optimality of the solution of the decomposed problem. In chapter 4, some 
specific information regarding the computer implementation of the EP design 
formulations is given, and results of using the two new formulations on test 
problems are obtained and compared with a commonly used method. The test 
problems for the EP design formulations range from a simple ten-bar truss to a 
high-speed civil transport having 348 design variables. The overall results are 
summarized, and suggestions for future work axe given in chapter 5. 



Chapter 2 


Equilibrium Programming Background 


In this chapter, the development of equilibrium programming and its previ- 
ous uses in engineering design are reviewed. The mathematical statement of 
an equilibrium programming problem, the necessary condition relations, and 
a solution existence theorem are stated. Some properties of an equilibrium 
programming solution and some solution methods are also given. 

2.1 Equilibrium Programming Background 

Equilibrium programming (EP), or non-cooperative game theory, was devel- 
oped in an operations research setting. The first proof of the existence of the 
solution to an equilibrium programming problem, called an equilibrium point, 
is due to Nash [34]. In Nash’s treatment, the subproblems assume the roles of 
players trying to maximize their individual pay-off functions in a game, and 
the domain for each subproblem is a fixed set of strategies that each player can 
utilize. The formulations of the subproblems are generalized in Debreu [10] so 
that the domains of each subproblem are functions of the design variables of the 
other sub problems. Debreu proves the existence of the equilibrium point, sub- 
ject to certain restrictions, assuming that the feasible domains are nonempty, 
closed, and bounded regions. Zangwill and Garcia [68] generalize the theorem 
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of Debreu, and prove existence of equilibrium points under weakened feasibility 
assumptions. 

Although there have been applications of equilibrium programming to 
economics, game theory, and network theory ([15], pp. 112 - 196), there have 
been few applications of equilibrium programming in engineering design. Its 
use in multicriteria engineering design is described by Vincent [65]. In mul- 
ticriteria design, tradeoffs between possibly conflicting criteria are required so 
that the design is acceptable to all the designers involved in the design process. 
Drawbacks in using a simplistic equilibrium programming approach to multi- 
criteria design (i.e., allowing each designer to independently optimize his own 
design objective function) are nonconvergence of the design, or convergence to 
a design that is inferior to other possible designs for all the designers. Im- 
position of an overall coordinating management of the design process may be 
required to obtain satisfactory designs. One coordination method for the game 
theory approach to multicriteria design is developed in [42] and is applied to 
illustrative structures-controls multicriteria design problems in [41]. Thus, an 
equilibrium programming formulation that includes some form of coordination 
appears necessary for equilibrium programming to be useful in engineering de- 
sign. In the EP design formulations derived in the present dissertation, the 
minimum-weight structural optimization problem, which is essentially a NLP 
problem, is reformulated as an EP problem. Coordination of the EP subprob- 
lems so that the solution is optimal can be important. 
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2.2 Features of Equilibrium Programming Problems 

The theory of equilibrium programming provides a framework to analyze multi- 
ple, interdependent nonlinear programming problems. Following Zangwill and 
Garcia [68], equilibrium programming is a generalization of nonlinear program- 
ming (NLP) which can be personified as having M decision makers (which 
may be implemented as search algorithms) that interact in a system. Each 
decision maker has a NLP subproblem to solve, and an independent set of 
design variables to control. The mathematical statement of this problem fol- 
lows. The design variables controlled by decision maker i are denoted as Xj, 
the design variables of all M decision makers are denoted as x = (xi, . . . , x A /), 
and all design variables not controlled by decision maker i are denoted as 
kj = (xi, . . . ,Xt_i,Xi+i, . . . ,x M ). Decision maker i has an objective function 
to minimize, /i(xj,Xj), while satisfying a set of constraints. Thus, the mathe- 
matical description of equilibrium programming is: 

min /i(Xi,Xi) 

subject to: g i (x i ,x i )<0 (2.1) 

hi (xj,Xj) = 0 

for the i — 1, . . . , M interacting NLP subproblems. The variables following the 
comma in any of the functions in statement (2.1) are treated as fixed parameters 
in that subproblem. Thus, in the NLP problem of decision maker i, the design 
variables from other decision makers x, enter as parameters that account for 
the coupling of the subproblems. 

Any value x that is a solution to all the NLP subproblems represented 
by statement (2.1) is called an equilibrium point. There may be numerous equi- 
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librium points of an equilibrium programming problem. In a manner similar 
to NLP, first-order necessary conditions are satisfied at an equilibrium point 
subject to a constraint qualification. Thus, at an equilibrium point, there exist 
Lagrange multipliers (A i,//*) such that the following conditions are satisfied: 

dfi{xi,Xi) , /x , T 3g i(xi,Xi) 

ax, + ax, 

. xT ^ h i(Xi,Xi) _ 

+ dXi 

gt (x»,Xi) < 0 (2.2) 

hi(xj, Xj) = 0 

A, > 0 

(A i ) T gi(x i ,x i ) = 0 

for i — 1, . . . , M. The conditions governing the existence of solutions (x, A, /x) 
that satisfy the necessary condition relations (2.2) and solve the equilibrium 
programming problem (2.1) are discussed in the next section. 

2.3 Existence of an Equilibrium Point 

The satisfaction of a constraint qualification is required for both the existence 
of a solution to the first-order conditions represented by statement (2.2) and 
for the existence of an equilibrium point. One form for the constraint qual- 
ification is given in [68]. The constraint qualification is satisfied if, for all x 
feasible to the NLP subproblems represented by statement (2.1), and for every 
subproblem i: 1) the vectors dh l j(x t ,x i )/dx l for all components j of h, are 
linearly independent, and 2) there is at least one solution z t to the relations: 


M Zi < o 

dxi 


(2.3) 


where g* is the vector of inequality constraints in subproblem i that axe active 
at x. With regard to the inequality constraints, this constraint qualification 
essentially states that it is always possible to move into the interior of the fea- 
sible region from a point on the boundary of that region. However, because 
constraint qualification relations (2.3) must be satisfied individually by each 
EP subproblem, these requirements on the constraint functions in EP are more 
restrictive than in NLP. For example, if the constraints from two EP subprob- 
lems are given by gi(-,x 2 ) and g 2 (x2, •), then relation (2.3) cannot be satisfied 
for subproblem i = 1 because {dg\/dxi) z x = 0. However, relation (2.3) may 
be satisfied for this example when the constraints and design variables of the 
subproblems are combined within a single NLP problem (i.e., (dg* /dx) z < 0 
where g* = (gf.gS) and x = (x!,x 2 )). 

A very general theorem for existence of an equilibrium point that is the 
solution of problem (2.1) is given in [68]. In this theorem, continuity, but not 
differentiability, of the objective and constraint functions is required. Other 
conditions for existence are: 1) the functions satisfy constraint qualification 
relations (2.3) (actually only a weakened form of the constraint qualification is 
required); 2) the feasible region is bounded with at least one feasible point x 7 for 
which g i (x',x i ) < 0 and h,(x 7 , x t ) = 0 for every feasible point x; 3) the func- 
tions fi(x it x t ) and g lJ (x z ,x- l ) are convex in x,; and 4) the functions /i tJ (xj, x t ) 
are linear in x, and, for a given i, have linearly independent gradients. The con- 
vexity and linearity restrictions on (/,, g t ) and h t , respectively, may be relaxed 
and a solution to the necessary condition relations (2.2) will still exist. How- 
ever, the solution may not be an equilibrium point. Other versions of the EP 
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existence theorem are reported in [68]. The sufficient conditions for existence 
of an equilibrium point are more restrictive than those required for existence 
of an optimal NLP solution. In addition to the more restrictive constraint 
qualification, a NLP solution exists with either condition 2) given above, or 
with condition 2) replaced by the coercive assumption defined by /(x) — + +oo 
as [x| — » oo ([69], p. 363). Some interesting equilibrium point properties are 
described in the next section. 


2.4 Some Equilibrium Point Properties 

Although the differences between an equilibrium point and an optimal point 
(i.e., the solution of a NLP problem) may appear slight, they are important. 
The following equilibrium point properties, summarized from [15], illustrate the 
differences. Several examples of these differences can be found in the reference. 

An equilibrium point is, in general, different from an optimal point, 
even if the same constraints are satisfied and each EP subproblem has the 
same objective function. This difference can occur because the coupling of the 
constraint derivatives in the respective necessary condition relations is gener- 
ally weaker for an EP formulation than for a NLP formulation. An example 
illustrating the differences follows. Suppose that there are two EP subproblems 
defined by the statements 


min 2x\ — x<i 

X\ 

subject to: 0 < X\ < 1 (2-4) 


—X\ + X2 < 0 
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and 


min 2xi — x 2 

subject to: 0 < x 2 < 1 (2.5) 

—X\ + x 2 < 0 

These subproblems have the same objective function and a common constraint. 
All the points on the line segment connecting (0,0) and (1,1) are equilibrium 
points. If any point on this line segment is attained during the solution process, 
the stability property of an equilibrium point ([15], p. 84) will prevent any 
movement from this point, even though there are neighboring solution points 
that would decrease the objective functions of both subproblems. Suppose 
that the constraint relations of statements (2.4) and (2.5) were combined into 
a single NLP using their common objective function, and with design variables 
(xi,x 2 ). Then the only optimal point for this NLP problem would be (0,0). 
EP formulations for modeling a system enable the M subproblems to have M 
different, and possibly conflicting, objective functions. Thus, in general, the 
equilibrium points of an EP formulation of a system will differ from the optimal 
points of an NLP formulation of the system that uses the same constraint 
relations but has only one objective function. 

Additional constraints can affect EP solutions differently than NLP so- 
lutions. In NLP, additional constraints generally increase the value of the 
objective function. However, in EP it is possible for additional constraints 
to force a coordination or cooperation of the subproblems that reduces the 
objective functions of all the EP subproblems. In the previous example, if 
the equilibrium point (0.75, 0.75) were found during the solution process, the 
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objective functions for both subproblems would be 0.75. If the additional con- 
straint X2 < 0.5 were added to subproblem (2.5), the objective functions at the 
new equilibrium would be reduced to a value no larger than 0.5. The stability 
property of an equilibrium point mentioned previously states that a solution 
to an EP problem does not change for a perturbation of the design variables 
of a single subproblem from equilibrium values. Under certain restrictive as- 
sumptions, the equilibrium point is unique, and therefore the stability of the 
equilibrium point is global [30]. If the equilibrium point is globally stable, the 
solution methodologies discussed subsequently are more likely to converge to 
the equilibrium point. 

2.5 Solution Methodologies 

Because the existence of an equilibrium point is independent of any particular 
solution method, an equilibrium point can be obtained in several ways. The 
most straightforward method is to solve all the subproblems sequentially in 
some predetermined cyclical order. When the solutions to all the subprob- 
lems do not change from the previous cycle, an equilibrium point has been 
reached. Although this method is used in the present dissertation, it may not 
converge, as demonstrated in [65]. Two approaches that improve the conver- 
gence characteristics of a sequential solution method are: 1) approaches that 
modify the objective functions so that they are strictly convex (such as the 
proximal point algorithm [44], or the penalty approach of [62]); and 2) move- 
limit-control methods. Incomplete convergence of intermediate subproblem 
solutions may also be used to improve computational efficiency. Another varia- 
tion of the sequential solution method is to solve the subproblems in a sequence 



21 


that is determined by information generated during the solution process. Some 
subproblems may even be omitted from a particular cycle. However, all the 
subproblems would need to be present in the last cycle to ensure convergence 
to an equilibrium point. 

In another solution method described in [15] (p. 97), the equilibrium 
point is found by converting the necessary condition relations (2.2) for each 
of the subproblems into a so-called Kuhn-Tucker equation set. The set of 
equations can then be solved using any nonlinear equation method. Homotopy 
methods are described in [15] as one method for solving these equations. 

In yet another solution method, applicable if parallel computation is 
available, the M subproblems can be executed on M processors in an asyn- 
chronous manner. Since the other subproblems communicate with subproblem 
i only through the parameters x, , the processors are effectively decoupled. Any 
update of x* during the solution process on processor i can be made immedi- 
ately available to the other processors. This method may only be applicable 
if the functions in the problem have favorable properties [30], and may require 
additional controls on the solution process to ensure convergence to an equilib- 
rium point. For example, relaxation techniques are used as a control method 
to improve convergence properties in [6]. 



Chapter 3 


Development of Equilibrium Programming 
Structural-Design Formulations 


In this chapter, the mathematical definition of a NLP-based, optimal struc- 
tural sizing problem is stated as a point of departure for the development of 
equilibrium programming structural-design formulations. Only a single load 
case is used in the following development; extension to multiple load cases is 
straightforward. Four EP structural-design formulations are then developed. 
The formal definition of the EP design variables, x*, will be given for each 
formulation. This formal definition may include both true design variables 
that can only be determined during the solution of the minimization problems, 
and behavior variables that can be determined from an analysis within the 
subproblem after the minimum is found. 

The first two EP design formulations developed herein were initially de- 
scribed by the author in [51]. Their primary purpose is to further acquaint 
the reader with equilibrium programming concepts, and to show that two com- 
monly utilized structural- sizing methods are, in actuality, EP formulations. 
Thus, EP formulations are currently being used to improve the computational 
efficiency of structural sizing. In these formulations, the EP subproblems con- 
sist of a structural-sizing subproblem (or subproblems for the first formulation), 
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and a structural-response subproblem for each load case. The first equilib- 
rium programming structural-design formulation developed herein considers 
only stress and local buckling constraints, and side constraints. This formula- 
tion is shown to be equivalent to the method of fully stressed design for rod 
and membrane elements. The second formulation to be described provides for 
optimal designs with displacement, stress and local buckling constraints, and 
side constraints. This formulation is shown to be equivalent to the NLP-based 
approach to optimal structural sizing using a first-order Taylor approximation 
of the structural response for a rapid analysis. 

In the third EP design formulation to be presented, an EP subproblem 
is developed that performs an approximate sensitivity analysis for a structural- 
response subproblem using the results of a single finite element analysis and the 
sensitivity derivatives of a previous iteration. Thus, the cost for this approxi- 
mate sensitivity analysis is negligible. This formulation was first reported by 
the author in [52]. Similar expressions for approximate sensitivity derivatives 
were developed using an ad hoc method in [31] for use in trajectory optimiza- 
tion, and the present formulation generalizes these expressions and provides a 
formal basis for their derivation. 

In the fourth EP design formulation developed herein, a novel substruc- 
turing technique is utilized to decompose the analysis and sensitivity calcu- 
lations for structural design. The substructural analysis portion of the for- 
mulation described herein was first reported by the author in [53]. In this 
formulation, the structure is divided into substructures, and each substruc- 
ture has its structural response and sensitivity derivatives determined by a 
structural-response subproblem. The structural sizing and the coordination 
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of the structural-response subproblems are determined by a single structural- 
sizing subproblem. 

3.1 Optimal Structural-Sizing Problem Statement 

The mathematical descriptions of finite-element-based, linear structural analy- 
sis, and optimal structural design are given in this section to provide a point of 
reference for the ensuing discussions. A single structural load case is assumed 
in the descriptions. 

A minimum potential energy formulation can be used in a finite-element- 
based structural analysis to calculate the structural displacements and stresses 
that are required to evaluate a design. Given the structural arrangement, the 
sizes for all the structural elements, a discretization of the structure into finite 
elements, and a set of external forces on the discretized structure, the correct 
structural displacements are those that minimize the potential energy of the 
structure. Thus, the structural response is the solution to the unconstrained 
NLP problem given by 

min (l/2u T Ku-F T u) (3.1) 

where the first term in (3.1) is the strain energy, and the second term is the work 
due to external forces. The vector u is a vector of nodal displacements, and F is 
a vector of external nodal forces (see the Nomenclature section for a discussion 
of notation and a full list of symbols). The domain of u is the entire space R 71 
(where n is the number of displacement degrees-of- freedom), but the boundary 
conditions on u are assumed to be incorporated in the global stiffness matrix K. 
The necessary condition relations for the unconstrained minimization problem 
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represented by statement (3.1) are simply the linear equations: 

Ku = F (3.2) 

Once the displacements u have been determined from either of these two state- 
ments, stresses a (or stress resultants N) within every element of the structure 
can be calculated from the displacements, the element strain-displacement re- 
lations, and the element constitutive relations. 

A nonlinear programming method that is used for optimal structural 
design can be personified as having one decision maker (which may be imple- 
mented as a search algorithm) with control of a set of design variables given by 
a vector v. The goal of the decision maker is to minimize an objective function 
W (v) while satisfying a set of constraints. The mathematical description of a 
NLP problem is: 

nun W(v) 

subject to: g(v) < 0 (3.3) 

h(v) = 0 

where g(v) are inequality constraints - which could include simple bounds on 
the design variables - and h(v) are equality constraints. A common choice for 
the objective function W(v) for structural sizing is the weight of the structure. 
The structural-sizing design variables v, referred to herein as sizing variables, 
can be the dimensions of the individual elements that explicitly contribute to 
weight, such as beam dimensions, skin thicknesses, and stiffener dimensions and 
spacing; or they can be variables which affect the weight in an indirect manner, 
such as the orientation of fibers in a composite structure. The constraint func- 
tions considered in the present study are side constraints (such as minimum 
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gauge) on the sizing variables g s (v), the local buckling and stress constraints 
g <T (v, <t(v)), and the displacement constraints g u (u(v)). The displacement 
constraint functions g u (u(v)) are assumed to have no direct dependence on 
the sizing variables v, but the functional form u(v) indicates an indirect de- 
pendence on v through the structural analysis. The stresses <x(v), which are 
shown to depend on v in the constraint functions, can have several forms. 
For example, one can write <x(v) = cr u (v, u(v)) = <x N (v, N(v, u(v))) to show 
that the functional form for stresses can depend directly on displacements u, 
or indirectly on displacements through the stress resultants N. Very often the 
side constraints, the displacement constraints, and the stress constraints reduce 
to simple bounds on the sizing variables, the displacements, and the stresses, 
respectively. Thus, the NLP approach to structural sizing that utilizes finite- 
element-based structural analysis can be summarized by: 

min W(v) 

subject to: g*(v) < 0 (3.4) 

g u (u(v)) < 0 

g CT (v,<r u (v,u(v))) < 0 

where the displacements u(v) are found from statement (3.1) or equation (3.2). 
The necessary conditions for the structural-sizing problem (3.4) are the same 
as given in (2.2) restricted to a single subproblem. 

3.2 Structural-Design Formulation 1 

In this section, a simplistic approach to defining an equilibrium programming 
formulation for structural design is described, its limitations are outlined, and 
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modifications that overcome the limitations are developed. In this approach, 
one EP subproblem is a structural-response subproblem that is defined by 
the minimization problem of statement (3.1). In this subproblem, denoted as 
subproblem 1, the design variables are the displacements (i.e., Xi = u), and 
the sizing variables v are treated as fixed parameters. By parallel reasoning, 
a second subproblem, denoted as subproblem 0, could then be identified as a 
structural-sizing subproblem represented by statement (3.4) in which the design 
variables are the sizing variables (i.e., x 0 = v), and the displacements u are 
treated as fixed parameters. The shortcoming of this EP formulation is that 
an equilibrium point may not exist because constraint qualification relations 
(2.3) cannot be satisfied for some constraints. For example, maximum stress 
constraints for rod elements having cross-sectional areas A as design variables 
(i.e., x 0 = v = A) are given by (x) = cr(u) — <x max < 0 . These constraints 
depend only on the displacements (xi), and not on the areas (xo); so there is no 
z 0 that will satisfy the inequality in relation (2.3). Thus, a more sophisticated 
formulation is required. 

A formulation that uses an alternate form of the stress constraints may 
enable satisfaction of the constraint qualification relations. Since functions for 
calculating the stress and buckling constraints can be constructed by using both 
the sizing variables and the stress resultants of a structure, a change of variable 
is made to utilize stress resultants in the constraint functions. Because the 
stress and buckling constraint functions depend explicitly on the sizing variables 
with this change of variable, the satisfaction of the constraint qualification is 
much more likely, but solution existence is still not guaranteed as will be shown 
in a subsequent example. 
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The stress resultants throughout the discretized structure are computed 
within the structural-response subproblems, and are symbolically represented 
by the equation 

N = N(v,u) (3.5) 

Thus, in this structural-design formulation, the structural-response sub- 
problem (subproblem 1) consists of a solution of the unconstrained minimiza- 
tion given by statement (3.1) (or its necessary conditions equation (3. 2)) for 
u followed by calculation of stress resultants N by equation (3.5). The design 
variables of the structural-response subproblem are defined as Xi = (u, N). The 
stress resultants N are used instead of displacements in formulating the stress 
and buckling constraints of the structural-sizing subproblem (subproblem 0). 
Thus, the structural-sizing subproblem in this equilibrium programming for- 
mulation is: 


min WYv) 
x 0 = v v 

subject to: v min — v < 0 (3.6) 

gSW^v.N)) < 0 

where the side constraints are shown as simple minimum gauge constraints in 
this subproblem. 

The method of fully stressed design [28] can be derived from this for- 
mulation if: 1) only one-dimensional rod and two-dimensional membrane finite 
elements are used; 2) one sizing variable is associated with each finite element 
having a stress constraint; 3) the stress constraints limit the maximum stress 
magnitude or von Mises stress; and 4) there are no buckling constraints. The 
structural-sizing subproblem can then be decomposed into a set of independent 
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elemental problems, one for each sizing variable and constrained element com- 
bination. The solution of these elemental problems is simple since the value 
of the sizing variable that makes a constraint active can be found analytically 
for each load case. For example, the elemental problem for the cross-sectional 
area of one-dimensional rod element j is solved by: 

Aj = max(A jVmm , \Nj\/a j<max ) (3.7) 

where max( ) chooses the maximum of its arguments, the structural-sizing 
design variable x 0j is defined to be the sizing variable Aj, and Nj is the axial 
force in the element j. Note that the quantities Nj are elements of the stress 
resultant vector N, and are also elements of xi and Xo- A solution method that 
alternates between solving the elemental sizing problems, and the structural- 
response subproblem leads to fully stressed design if it converges. 

As stated previously, the change of variables that recasts the stress and 
buckling constraints in terms of sizing variables and stress resultants makes the 
existence of an equilibrium point more likely, but not guaranteed. A simple 
example makes this statement clear. Assume a rod is fixed between two rigid 
walls and its temperature is increased; the design problem is to size the rod 
cross-sectional area A to minimize the weight and to satisfy a maximum stress 
constraint. The temperature change induces a strain which can be expressed 
as an equivalent external load that is a function of the stiffness. Thus, the 
equivalent external load (which is also the rod stress resultant N ) is calculated 
by the equality constraint hi = N - k AT A = 0 in the structural-response 
subproblem where the sign of N is defined to be positive for compression, k is 
a constant of proportionality, and AT is the temperature rise. The minimum 
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gauge and the stress constraints, which are calculated in the structural-sizing 
subproblem, are given by g 0 , 1 = A min - A < 0 and g 0 , 2 = N - Aa max < 0, 
respectively. These constraints are shown in figure 3.1. Since a minimum 
weight design corresponds to a minimum A that satisfies the constraints, two 
conditions can be identified in the figure. For a relatively low value of ATi, 
the fully stressed design approach given here will converge to x* shown in the 
figure. If AT 2 is above a limiting value, the constraints allow for no feasible 
region in the design space, a prerequisite for solution existence, and the fully 
stressed design algorithm would diverge. This example is severe because even 
a more sophisticated design method would fail for a large enough AT because 
there would be no feasible region in design space. 

Because structural design formulation 1 is a form of fully stressed de- 
sign, it shares the advantages and disadvantages of fully stressed design. The 
primary advantage is the simple nature of a structural-sizing subproblem that 
requires no derivatives and is easily decomposed into a set of independent ele- 
mental problems. The disadvantages are: 1) there is no mechanism to ensure 
satisfaction of the necessary conditions of the optimal structural-sizing problem 
(3.4) so that the resulting equilibrium point may not be an optimal point; and 
2) constraints, such as displacement constraints, that have no explicit depen- 
dence on the sizing variables are not considered. An EP formulation is desired 
which satisfies the optimality necessary conditions at an equilibrium point, and 
can satisfy displacement constraints. 
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Figure 3.1: Design space for a rod constrained between two rigid walls under- 
going a temperature change. 

3.3 Structural-Design Formulation 2 

If the EP subproblems are identified using the simplistic approach described at 
the beginning of the previous section, the constraint qualification relations (2.3) 
cannot be satisfied for any displacement constraint because the displacements 
would be fixed parameters within the structural-sizing subproblem. In addition, 
the constraint qualification relations may not be satisfied for certain stress and 
buckling constraints. Substituting approximate models, that depend explicitly 
on the sizing variables, for these displacements within the structural-sizing 
subproblem can overcome these difficulties. 

In EP structural-design formulation 2, the displacements u, which are 
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parameters independent of the sizing variables in structural-sizing constraint 
functions g(v, u), are replaced with a first-order Taylor series approximation 
given by 

u 4 (v,X!) = ^ (v- vj) + u (3.8) 

In equation (3.8), the matrix <9u/<9v may be viewed as a matrix of 
optimal sensitivity derivatives with respect to parameters v [55] of the dis- 
placements determined in the structural-response subproblem (subproblem 1). 
The value v x is the value of v utilized in subproblem 1 when the sensitivity 
derivatives are calculated. All the design variables for the structural-response 
subproblem are utilized in equation (3.8) since Xj = (u, d\i/dv, v x ). This ap- 
proximation satisfies the following two properties. First, the approximation 
depends explicitly on v so that constraint qualification relations (2.3) neces- 
sary for equilibrium point existence can be satisfied. Second, the approximation 
satisfies the conditions: u A = u and du A jd\ = du/dv at the equilibrium point 
x = x* where v = v x = v*. This second set of conditions ensures the opti- 
mality of the design given by the equilibrium point because the EP necessary 
conditions are then the same as the NLP necessary conditions for problem (3.4). 

Using the definition of equation (3.8), the structural-sizing subproblem 
is given by the following statement: 

min W (v) 

Xo = v 

subject to: v l < v < v“ 


g» < 0 

g u (u A (v,x x )) < 0 
g^Vj^v.uV.xO)) < 0 


(3.9) 
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where move limits v ; and v“ that are adjusted during the solution process have 
been introduced to ensure convergence. The structural-response subproblem is 
represented by unconstrained minimization problems given by statement (3.1) 
(or the necessary conditions given by equation (3.2)), and the following equa- 
tions (which can be formally treated as equality constraints) that determine 
the behavior variables du/dv and vy 

= [0] (3.10) 

= V 




Vi 


The first of the two equations in statement (3.10) is a sensitivity analy- 
sis equation using a direct method of sensitivity analysis ([19], p. 264). In this 
equation, the sizing variables v are assumed known, the vector u in the second 
term is treated as a constant, and contributions from the nodal forces F are ne- 
glected since F is assumed to be independent of v. In practice, the second term 
of this equation is often determined by finite differences. The second equation 
in (3.10) is a trivial identity which preserves the values of the structural-sizing 
variables at which the sensitivity analysis is performed so they can be exported 
to the other subproblems. 

EP structural-design formulation 2 is equivalent to a NLP approach to 
structural design with a first-order Taylor series for an approximate analysis 
(i.e., the method of [45] except that reciprocal variables are not used) if the 
equilibrium point is found by the following sequential steps: 1) solve the sub- 
problems represented by statement (3.1) (or necessary conditions represented 
by statement (3.2)) for u; 2) solve equations (3.10) for du/dv and vp, 3) use 
the quantities found in steps 1 and 2 in equation (3.8) to form approximate 



models for displacements u' 4 (v,x i ); 4) utilize the approximate models of step 
3 in the structural-sizing subproblem represented by statement (3.9) to solve 
for v; and 5) repeat steps 1 through 4 in a cyclic manner, with an algorithm 
defined to update the move limits, until the changes in the solutions from con- 
secutive cycles converge. Although first-order Taylor series are used in this 
development, other approximate models for displacements that depend on first 
derivatives could be used in place of (3.8) ( [19], pp. 211-219). 

3.4 Structural-Design Formulation 3 

The two EP design formulations described previously represent current prac- 
tices in optimal structural sizing. To further improve computational efficiency, 
a formulation is derived herein that computes approximate, updated sensitivity 
derivatives in an EP subproblem. The cost of computing these approximate 
derivatives is essentially only the cost of the structural analysis, and the eval- 
uation of the design constraints. 

The approximate sensitivity analysis is derived as a correction, or up- 
date , to previous values for the sensitivity derivatives. The sensitivity-derivative- 
update subproblem is described for a general scalar function /(v), where v is 
a vector quantity, as follows. The values of / and the gradient of /, df/dv, 
are assumed to be known at a previous value of v. These quantities are de- 
noted as / 0 id, df 0 \d/dv, and v o]d , respectively. The value of / at the new value, 
v = v new , is also known and an approximation to the gradient of / at v = v new 
is sought; this approximate gradient is denoted as df new /dv. The difference be- 
tween <9/new/d v and df 0 \d/dv is given by the gradient updating vector a (also 
called an update vector, with components that axe called updates herein), and 
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the relation among these quantities is expressed by the equation 


5/new 5/old , 

-0T- = ^T + a 


(3.11) 


Using a criterion similar to that utilized in [11] (p. 171) in calculating 
least-change secant updates of the Jacobian matrix for solving simultaneous 
equations, the update vector a is chosen to be the vector that has the small- 
est (weighted) magnitude, and that also satisfies a second-order Taylor series 
relation between the quantities at v = v new and v = v 0 id- The constrained 
minimization problem that expresses these conditions is: 


min T(a k w k f 

a/c , 
k 

subject to: /new = /old + + ^ a k)( v new,k - u old,fc) (3.12) 

where w k are user-defined weighting parameters. The equality constraint in the 

minimization problem given by statement (3.12) indicates that the quantities a k 

are approximations to X) 5 2 //5u*5u/ (v new i — v 0 ] d j). Thus, the accuracy of the 

i 

sensitivity derivatives determined by the update vector is expected to degrade 

as the quantity |v new - v o]d | increases in value. Letting Av k = ^new ,k ^o\d,k 

and A / = / new - f 0 \ d - E df o]d /dv k Av h , the solution to the minimization 

k 

problem given by statement (3. 12) can be found analytically using a Lagrange 
multiplier method 


Q, k — 2 


A/ Av k 


w k E A vf/wf 

l 


(3.13) 


A simplified version of this expression has been developed by an ad hoc 
method in [31] with application to the approximation of gradients for trajec- 
tory optimization. However, the update equations described herein are derived 
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independently using a formal methodology that allows for the general inclu- 
sion of the weights w k that are not considered in [31]. For identical weights 
w k , the update given by equation (3.13) is twice the Broyden update equation 
described in [11] (p. 170), and is the same as that derived in [31]. The advan- 
tage of including the general weighting terms is described subsequently. In the 
present study, the function / represents a component of the constraint function 
vector g. 


3.4.1 Structural-Design Formulation Using Updated Sensitivity Deriva- 
tives 

The method of approximate, updated sensitivity derivatives is denoted as EP 
structural-design formulation 3 (or simply formulation 3) herein. This for- 
mulation is a modification of design formulation 2 described previously in 
which update subproblems (3.12) are utilized during the solution process in 
place of computing sensitivity derivatives such as (3.10). In subproblem 0 (the 
structural-sizing subproblem), the design variables are Xo = (v, v 0 , go, <3go/<Tv). 

All but the first of which are actually behavior variables in that subproblem. 
Formally, this subproblem can be represented by the following statement: 


nun W(v) 

subject to: ^-(v - Vi) + gi 

(JV 

v l < 

g s (v) 

Vo 

go 


< o 

v <v u (3.14) 

< 0 

= Vj 

= gi 



where the general constraints are linearized, and move limits v l and v u that 
are adjusted during the solution process have been introduced to ensure con- 
vergence. The last three equalities in problem (3.14) are required to define 
data items that are saved to be utilized in the structural-response subproblem. 
Subproblem 1 (the structural-response subproblem) has design variables given 
by xi = (u, gi, vi, dgi/dv), and most of these quantities can be computed after 
solving for the first member of Xi (i.e., u). The subproblem statement is given 
as: 

min ^l/2u r K(v)u — F T u^ 

subject to: gi = g(v,u) (3.15) 

Vi = v 
SEN S (choice) = 0 


where the expression SEN S {choice) = 0 represents the conditions: 


K(v) 


du 

dv 


dg ij 

dvi 


dK(v)u 

&v 

gj ((v + Auie7), (u + du/dv A^e^)) - gij 

Avi 


(3.16) 


when “exact” sensitivity analysis is desired, and 

dg u dgoj . nf , ( dgoA 2 9ij - 9oj ~ dgojdv (v - v 0 ) 

~dA~dv t Z[Vl Vo ’ l} V dvi J Zi(vi - votndgoj/dvtf 


(3.17) 


when approximate sensitivity analysis is desired. Several implementation fea- 
tures of this formulation should be noted. The potential energy" minimization 
in (3.15) is a formal statement, and may replaced with the necessary conditions 
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of equation (3.2). The equality constraints in (3.15) are formal definitions to 
define quantities needed in the structural-response subproblem. As in the de- 
termination of the move limits, an additional algorithm is required to determine 
whether statement (3.16) or statement (3.17) is utilized for SEN S(choice) = 0 
in statement (3.15) at each stage in the solution process. Because the update 
minimization subproblems have analytic solutions given by equation (3.13), 
these solutions are incorporated as constraints in the structural-response sub- 
problem using equation (3.17) instead of specifying them as independent sub- 
problems. Finally, the weights in the update subproblems have been specified 
in equation (3.17) as the inverse of the gradients, Wi — dgo/dvi. Thus, the 
updates to the sensitivity derivatives determined from the subproblems given 
by statement (3.12) minimize the relative changes of the sensitivity derivatives 
from their previous values, and the update to each component of a gradient 
vector is related to the component magnitude. One advantage of this weight- 
ing method is that the sparsity pattern of sensitivity derivatives is maintained 
since an entry that is zero or small in magnitude remains zero or small in mag- 
nitude after updating. This method of weighting was also found to provide 
better convergence than using equal weights in some preliminary optimization 
studies. More details concerning the computer implementation of this design 
formulation are discussed in the next chapter. 

3.5 Design Formulation 4 

In this section, an EP design formulation is described that utilizes decom- 
position methods when computing the structural response and the sensitivity 
derivatives. In this formulation, the structure is divided into substructures, and 
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each substructure has its structural response described by a structural-response 
subproblem. The structural sizing and the coordination of the structural- 
response subproblems are determined by a single structural-sizing subproblem. 

Substructure-based design methods have been utilized in [25] and [26], 
but the analyses in these references were not based on substructure princi- 
ples. Traditional (i.e., super element), substructure-based sensitivity methods 
for structural optimization have been used to form a reduced basis for approxi- 
mate analysis in [36] , have been derived using the adjoint sensitivity method in 
[1], and have been combined with kinematic constraints to reduce the number 
of interface degrees of freedom in [37]. The adjoint method of substructure 
sensitivity analysis of [1] is generalized to multiple levels of substructures in 
[35]. An alternate method of using substructures for performing structural 
analysis is developed in [12]. The method of reference [12] is essentially a 
substructure-based, hierarchical decomposition method in which the response 
of each substructure is determined by independent subproblems. Coupling con- 
straints that determine the interface forces required to enforce the compatibility 
of common interface nodes form a coordination problem. Decomposition into 
substructures that are not restrained from rigid-body motions yields subprob- 
lems that are not strictly convex. Nonconvexity causes difficulties, as noted 
in [17], that can be handled using special techniques for solving both the sub- 
problems for each substructure and the coordination problem. A modification 
of this alternate substructure analysis method, described in [13], incorporates 
penalty functions to ensure strict convexity of the subproblem objective func- 
tions. 


The EP design formulation of the present section is similar to the 
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method of [12], but strict convexity is not required in this approach. In the 
present derivation, each design constraint is assumed to be local so that it 
depends on the sizing variables and the structural response of a single sub- 
structure. Thus, stress constraints, local buckling constraints, and constraints 
that axe functions of displacements within single substructures may be con- 
sidered in this approach. In the next subsection, the decomposition of the 
structural analysis is discussed utilizing a two-substructure example. Then 
the derivation is extended to an arbitrary number of substructures, and the 
substructure-based EP design formulation is developed. 

3.5.1 Derivation of the Substructure Analysis Method 

The method of determining the structural response is presented for a struc- 
ture that is decomposed into substructures, each of which is assumed to have 
a linear elastic response with no zero-strain-energy motions possible except for 
rigid-body motions. To permit a simple exposition of the method, it is initially 
derived using two substructures. Four salient features characterize the present 
method for structural analysis using substructures. Firstly, when the substruc- 
ture is not fully restrained from rigid-body motions, the structural response 
is decomposed into the rigid-body motions (referred to as “modes” herein) of 
the substructure, and displacements orthogonal to the rigid-body motions (i.e., 
the elastic deformations). Secondly, an augmented stiffness matrix is formed 
for each of these substructures. These stiffness matrices are symmetric, and 
can be factored independently and, computationally, in parallel. Thirdly, a 
structural-response coordination problem determines the forces between the 
substructures, and the magnitude of the rigid-body modes. The structural- 
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response coordination problem requires the factored, stiffness matrices of the 
substructures, and results in a system of linear equations with its order equal to 
the number of shared degrees-of- freedom between substructures plus the total 
number of substructure rigid-body modes. Lastly, once the forces between the 
substructures are determined, the displacements orthogonal to the rigid-body 
modes can be determined. 

The simple wing structure finite element model in figure 3.2 is utilized to 
demonstrate the derivation of the analysis decomposition method. The model 
shown in the figure is decomposed into two substructures, the second of which 
is unrestrained and thus has six rigid-body modes. The linear elastic structural 
response for the entire wing is determined using a minimum potential energy 
formulation. This response is the solution to the unconstrained minimization 
problem given by statement (3.1) having the necessary conditions given by 
equation (3.2). In this section, the nodal displacements are denoted using the 
lower case (i.e., u). When the displacements of a substructure can be decom- 
posed into elastic deformations and rigid-body modes, the elastic deformations 
are denoted by the upper case (i.e., U). Thus, the displacements of the entire 
wing are denoted by the nodal displacement vector u, and the displacements of 
the two substructures are given by the vectors u x and u 2 . The external loading 
for the entire wing is given by the nodal force vector F, and the decomposition 
of the external loading is given by vectors Fi and F 2 for the two substructures. 
The dimensions of each of these vectors is the number of nodal degrees of free- 
dom of the appropriate substructure. The nodal displacements at the interface 
of the two substructures must be equal for the two substructures to be compat- 
ible. Those compatible displacements at the interfaces of the two substructures 
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Figure 3.2: Example wing structural finite element model decomposed into two 
substructures. 

are specified in a predetermined order with signed Boolean matrices denoted 
by Bj and B 2 , where the nonzero entries of Bj are equal to +1 and the nonzero 
entries of B 2 are equal to —1. The constraints of compatible displacements at 
the interface between substructures is then expressed by B 1 Ui+B 2 u 2 = 0. The 
decomposition of the externally applied nodal forces at the interface nodes is 
arbitrary as long as the sum of the forces applied to the interface nodes of each 
substructure is equal to the actual externally applied forces at these nodes. 

In a substructural decomposition, a stiffness matrix can be formulated 
for each substructure. These stiffness matrices are relative to the substructure 
nodal displacement vectors u, . and are denoted by Kj. Thus, a constrained 
minimization problem that is equivalent to the problem given by statement 
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(3.1), and uses the substructure nodal displacement vectors is 

min 

(ui,u 2 ) l 

-Ff Ul -F*u 2 ) (3.18) 

subject to: BiUi + B 2 u 2 = 0 

The constrained minimization problem given by statement (3.18) has a separa- 
ble objective function with coupling constraints, and appears to be a candidate 
for a hierarchical decomposition into two distinct minimization problems hav- 
ing design variables Ui and u 2 , respectively. The lower-level subproblems for 
the displacements using such a decomposition are 

min ^ ufKxUi - Ffuj + A t BiUi (3.19) 

Uj 2 

and 

min ^ulK 2 u 2 — Flu 2 -1- A T B 2 u 2 (3.20) 

u 2 2 

where the goal coordination approach described in [32] (p. 240), and sum- 
marized in the first chapter of the present dissertation, is utilized with the 
common coordination inputs A. To find the minimum of problem (3.18) using 
goal coordination, the coordination inputs would need to be determined by a 
coordination problem that is typically the dual formulation of statement (3.18). 
However, the minimization problem for (3.20) is insoluble due to the rank de- 
ficiency of matrix K 2 . An EP explanation for the insolubility of subproblem 
(3.20) can be developed based on the conditions for the existence of an EP 
solution. One of the sufficient conditions for the existence of an equilibrium 
point requires a bounded solution space for the subproblems. However u 2 is 
not bounded as seen from the following argument. Because substructure 2 has 
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six rigid-body modes, a rigid-body matrix R 2 can be formed that consists of 
the six rigid-body modes as column vectors. This rigid-body matrix satisfies 
the condition K 2 R 2 = [0], thus the vector R 2 a 2 added to any displacement u 2 
will not change the substructure potential energy if the vector ot 2 satisfies the 
scalar equation F 2 R 2 a 2 = 0. Since the vector a 2 is not otherwise restricted, 
it (and therefore u 2 ) is unbounded. Thus, the sufficient conditions for solution 
existence are not satisfied. To overcome the difficulties due to rank deficiency 
of substructures having rigid-body modes, a pseudo inverse of the matrix K 2 
is used in [12] when solving for u 2 in the necessary conditions for subproblem 
(3.20). The coordination inputs A, and the magnitude of the rigid-body modes 
are determined in a coordination subproblem. Alternate ways of modifying 
subproblem (3.20) to make its solution bounded are to add terms to make the 
objective function strictly convex using the proximal point algorithm as in [7], 
or using the separable penalty method of [62]. Another penalty method that 
can convexify the problem is reported in [13]. 

From the preceding discussion, any decomposition of statement (3.1) 
should explicitly account for the rigid-body modes of the substructures. One 
such decomposition method follows. Let u 2 = U 2 +R 2 o: 2 where U 2 satisfies the 
supplementary condition R 2 U 2 = 0 (i.e., U 2 is orthogonal to the rigid body 
modes). The minimization problem given by statement (3.18) then becomes 

, min % (iufK lUl + W 2 U 2 
(u!,U 2 ,a 2 ) z 

-F[ Ul - F 2 (U 2 + R 2 a 2 )) (3.21) 

subject to: BxUi + B 2 (U 2 + R 2 a 2 ) = 0 

r!u 2 - 


0 



Using a hierarchical decomposition with goal coordination, three lower-level 
subproblems can be formulated. The first is identical to statement (3.19), the 
second is 
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min (^U[K 2 U 2 - F[U 2 + A T B 2 U 2 ) (3.22) 

U2 / 

subject to: RjU 2 = 0 

and the third is 

min — F^R 2 a 2 + A T B 2 R 2 a 2 (3.23) 


The objective function of subproblem (3.22) is strictly convex relative to 
the subspace (U 2 | R^U 2 = 0}, and thus the subproblem will have a bounded 
solution. However, subproblem (3.23) will have no finite solution unless A sat- 
isfies the restriction — F^R 2 -I- A T B 2 R 2 = 0 r . This restriction on A becomes 
a constraint in the standard goal coordination subproblem that maximizes the 
Lagrangian function (as a function of A) in the dual formulation of statement 
(3.21). This restriction is also used in the conjugate projected gradient ap- 
proach to solving the coordination subproblem in [12]. One expression for the 
coordination subproblem can be found from the following necessary conditions 
for the structural-response problem given by statement (3.21) 


K1U1 — F 1 + Bf A = 0 


(3.24) 


U 2 ] 


F 2 - B^A 



0 


R2 (F 2 — 63 A) = 0 


(3.25) 

(3.26) 


B1U1 + B 2 (U 2 + R 2 q: 2 ) — 0 


(3.27) 


where A are the Lagrange multipliers for the first constraint in statement (3.21), 
fj, are the Lagrange multipliers for the second constraint in statement (3.21), 
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and the matrix K 2 is the augmented stiffness matrix given by 


K 


2 


k 2 

r 2 1 


[0] J 


(3.28) 


Because the matrices B, in necessary conditions (3.24-3.27) are Boolean, the 
Lagrange multipliers A are simply the interface forces between the substruc- 
tures. Using equation (3.25), the condition RjK 2 = [0], and the fact that 
the columns of R 2 are linearly independent, equation (3.26) implies that the 
Lagrange multipliers (i are equal to zero. Note that equation (3.24) is the 
necessary condition equation for subproblem (3.19), and equation (3.25) is the 
necessary condition equation for subproblem (3.22). 


The structural-response coordination subproblem that determines A and 
a .2 is obtained by substituting the solutions of equations (3.24) and (3.25) into 
equations (3.26) and (3.27). This substitution yields the equation 


M 


A 


' BiK^F^ fB 2 |[0]lK 2 1 [Fi’|0 7 ’] T ' 

ct 2 


- R 2?2 


(3.29) 


where the matrix M is given by 
M = 


B 1 Kr 1 B[ + 

B 2 1 [0] 

KJ 1 

b 2 1 (0] 

T 

— B 2 R 2 

[ — (B 2 R 2 ) r 

[0] J 


(3.30) 


Note that the inverse of K 2 used in equations (3.29) and (3.30) exists because of 
the inclusion of the constraint that U 2 is orthogonal to the rigid-body modes. 
Thus, one approach to the calculation of the structural response (i.e., the nodal 
displacements) is the following steps: 1) factor the matrices K x and K 2 ; 2) use 
these factored matrices to formulate the matrix M and the right hand side 
of equation (3.29); 3) solve equation (3.29) for A and a 2 ; and 4) and solve 
subproblem (3.19) (or equation (3.24)) and subproblem (3.22) (or equation 



(3.25)) for Ui and U 2 , respectively. The explicit computation of matrix M is 
not necessary in the combined analysis and optimization approach described 
in the next section. 


47 


3.5.2 Derivation of EP Design Formulation 4 

In this subsection, the structural-optimization problem is decomposed into sub- 
problems that perform the structural analyses and sensitivity analyses of each 
substructure and a single subproblem that performs the coordination of the 
structural analyses as well as the optimization of the design. This formula- 
tion is denoted as EP design formulation 4. The decomposition method for 
structural analysis derived in the previous section is extended to multiple sub- 
structures. Then the design formulation is derived starting with a simultaneous 
analysis and design formulation. The simultaneous analysis and design method 
was initially investigated in [46]. It is utilized in [22] for linear structural anal- 
ysis, and in [20] for nonlinear structural analysis. In the present simultaneous 
analysis and design formulation, in which minimum weight is the design goal, 
both the structural displacements and the sizing variables are utilized as design 
variables, and the optimization constraints are the equations governing struc- 
tural response (treated as equality constraints) as well as the usual inequality 
constraints that ensure that the design meets the strength, buckling, and other 
design requirements. Formulation 4 has the following features: 1) the analysis 
and sensitivity derivatives of each substructure are independent and can be 
performed in parallel (although the overall optimization procedure is still iter- 
ative); 2) the resulting design is optimal; but 3) the structural response may 
not be compatible between the substructures until the design converges. 
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Structural Analysis using Multiple Substructures 

The substructural analysis method developed previously is generalized to mul- 
tiple substructures in this subsection. With multiple substructures, the use of 
precise, although somewhat cumbersome, notation is necessary to avoid confu- 
sion. Assume that there are M substructures, and that there are m substruc- 
tures (0 < m < M ) unrestrained from rigid-body motions if separated from the 
rest of the structure. The rii displacement degrees-of-freedom of substructure 
i are denoted by u, = (u h i, . . . , u itni ) T . The substructures unrestrained from 
rigid-body motions are ordered to be the last m substructures. These substruc- 
tures have displacements denoted by u, = U* + R,a, where R, is the matrix 
containing the rigid-body modes of substructure i, and Uj are displacements 
orthogonal to these rigid-body modes. The orthogonality relation for Uj is ex- 
pressed by RfUj = 0. There are n 1 independent relations governing the com- 
patibility of the substructures, and the resultant equivalencing of degrees-of- 
freedom at common interface nodes in the different substructures is represented 
by the set A of n 1 4-tuples defined by A = { (j, k,p , q) \ j < k and u hV = Uk, q }. 
In the definition of set A, redundancies in the equivalencing of degrees-of- 
freedom are omitted. For example, a degree-of-freedom shared by three sub- 
structures need only be equivalenced between two pairs of substructures, not 
all three pair-wise combinations. The compatibility constraint equations (i.e., 
the equations that enforce compatibility between the substructures) are then 
defined in the following manner. The r th compatibility constraint equation will 
depend on a degree-of-freedom in substructure j if either the first or second 
element of the r th 4-tuple in A is j. These compatibility constraint equa- 
tions are expressed explicitly by defining the signed Boolean matrices B j for 
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j — 1, . . . , M which have dimensions n 1 x rij. Matrix Bj has a 1 at location 
(r,p) if the r th 4-tuple in A has the components (j, ■ ,p , •), and it has a —1 at 
location (r, q) if the r th 4-tuple in A has the components (-,j, -,q). Otherwise, 
the entries in matrix B 7 are zero. Thus, the compatibility constraint equa- 
tions for all the substructures are symbolically represented by the system of n 1 
equations: 

M 

53 BjUj = 0 (3.31) 

i - 1 

The definition of the Boolean matrices B_, in the present section reduces to the 
previous definition of the Boolean matrices for the case of only two substruc- 
tures. 


Utilizing these definitions, the minimization problem for the structural 
response given by statement (3.21) generalizes to: 


M—m s I 

min (-ujKiUi 

(Ui, . . . ,UA/-mj i=l 
u M — ra+1? • * • i Uju, 

&M— rn+lj • • • i &m) 



subject to: 


+ E Ur KiU, - Ff (U, + iw) 

M—m M 

53 BiUi + 53 Bj(Uj + RjOLj) = 0 

z=l j=M—m + 1 

RjUj = 0 for i = M — m + 1, . . . , M 


(3.32) 


The necessary conditions for the minimization problem given by statement (3.32) 
are 


KjUj - Fi + Bf A = 0 


Ui 1 


Fi - Bf A 1 

. /*» . 


0 


Rf (Fi ~ BjX) = 0 


for i = 1, . . . , M — m 

for i = M — m + (3.33) 


for i = M — m + 1, . . . , M 
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along with the first constraint equation given in statement (3.32). Here A are 
the Lagrange multipliers corresponding to the first constraint equation in state- 
ment (3.32), fj, t are the Lagrange multipliers for the last constraint equations 
in statement (3.32), and the matrices K* for i = M — m + 1, . . . , M are the 
augmented stiffness matrices given by 



The second equation in necessary conditions (3.33) combined with the condi- 
tions RfKj = [0] for i = M — m + 1, . . . , M, and the fact that the columns of 
Ri axe linearly independent implies that the last equation in necessary condi- 
tions (3.33) can be replaced with the condition that the Lagrange multipliers 
/r t are equal to zero. The structural responses u, for i = 1, . . . , M — m, and 
(Ui,^) for i — M — m + 1, . . . ,M are found by solving the first and sec- 
ond equations of necessary conditions (3.33), respectively, after the vector A is 
determined. 

The vector A may be determined simultaneously with the rigid-body 
displacements a l for i = M — m + 1, . . . , M by a coordination subproblem. A 
coordination equation may be found by substituting the symbolic solutions for 
u i and Uj into the compatibility constraint equation in statement (3.32), and 
solving this equation simultaneously with the last equation in necessary condi- 
tions (3.33). The following matrix equation gives the coordination subproblem: 
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\ 


■ BiKr'F. + ES „- m+1 [B. | [0]]K-> [Ff I O^f ' 


&M—m + 1 


— m + 1 ^ M ~ m + 1 

M 

OCM-m+2 

= 



Ot-M 


1 

£ 

” p? 
i 

i 


(3.35) 


where the matrix M is given by 


E":”>B.Kr 1 Bf+ 

E"„- m+1 fB.|[0]lKr'[B.I[0]] T 



1 

§ 

§ 

CD 

1 

(Bm 

[0] 


[0] 

M —m+ 2 ) 

[0] 


[0] 




l 

— (BmR-m) t 

[0] 


[0] J 


l- ' ' I «■ ■» l l *■ J -J 

(3.36) 

Statements (3.35) and (3.36) are the generalizations of statements (3.29) and 
(3.30) for multiple substructures. These statements are given for completeness, 
but are not utilized in the following design formulation. 


Decomposition of the Structural- Design Problem 

The starting point for the design decomposition is a simultaneous analysis and 
design formulation. The constraints are assumed to have the form gi(vj, Uj) 
for i = 1, . . . , M in this derivation. The simultaneous analysis and design 
formulation is given by 

M 

min ^W^Vj) 

(vi,.. . ,v M , i=l 

- j HA/— rm 

U M — m+l> * • • U M , 

OCM-m+U • • • ) a Mi 
P'M — 771+1 ! t Pm) 

subject to: gj(vj,Ui) < 0 for i = l,...,M — m 



52 



for i = M — 

for i = 1 , . . . , M — m 

for i = M — m + 1, . . . , M 

(3.37) 

for i = M — m + 1, . . . , M 


This expression is simplified by the following steps: 1) specify that u* for 
i — 1, . . . , M— m and U* for i = M — ra+1, . . . , M (which are implicit functions 
of Vj and A) are known from solutions of structural-response subproblems (such 
as problems (3.19) and (3.22)), and thus identically satisfy the first two equality 
constraints in statement (3.37); and 2) for the remaining functions in (3.37), 
replace g,, u and U, with their first-order Taylor series approximations about 
the point (v*, a i} A). The design variables for the structural-sizing subproblem 
are given by x 0 = (v 1} . . . , v M , ctM-m+i > • ■ - , &m, A). These simplifications 
reduce statement (3.37) to the following structural-sizing subproblem 

M 

min Wi( v i) 

x 0 = (Vi,. ..,V M , t=l 

i • • • 5 Q M 7 A) 


subject to: 

v{ < Vj < v“ 

for 

i = 1,.. 

.,M 


gt + 

dvi x 

l) <9 A 

(A - Aj) < 0 

V< 

for 

i = 1, • • 

. ,M — m 


K .+ a&i 
dvi 

(Vi- 

A ,a, 

v°)+^i 
l) v , 

(A - \) 

,a t 







+ |^ (oci — Q!j ) < 0 

8a • v„A 

for 

i = M- 

-m + 1,.. 

■ > M 

M-m 

J2 Bi (uj + 

i= 1 

dn i . 

■X— (Vj 
dvi x 

v?) + 9u ' 

l> d\ 

(A - A.)) 

Vi 




(3.38) 
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M 


+ E + g 

au/ 


(Vi-vj) 


+ 


dX 


(A — A j) 3“ RjCkj^ — 0 
Rf(Fi - Bf A) = 0 


for i — M — m-\-\,...,M 


where move limits v- and v“ that are adjusted during the solution process 
have been introduced to ensure convergence. In this subproblem, the design 
variables A are coupling variables, and the compatibility constraints (the next 
to last constraint shown in statement (3.38)) axe coupling constraints. Be- 
cause of the relatively small number of these coupling variables and coupling 
constraints, a further decomposition of this structural-sizing subproblem into 
smaller subproblems could be attempted using the method of Ritter ( [29], pp. 
276-283), or the method of Ha [18]. However, no further decomposition of the 
structural-sizing subproblem is undertaken in the present dissertation. 

The structural-response subproblems, one for each substructure, deter- 
mine the structural displacements, the design constraints, and their derivatives 
with respect to the elements of x 0 . The design variables of the structural- 
response subproblems for substructures i = 1, . . . , M — m are x, = (u*, v°, A,, 
du l /dv l \^ , chij/dA| v . , g,, Sgi/dv^ , <9g,/dA| Vi ). The necessary condition 
form of these subproblems is 

KjU, = F, - Bf A 

v° = Vi 

Ai = A 

g i = gi(Vi,Ui) 
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for i — 1, . . . , M — m. The design variables of the structural-response subprob- 
lems for substructures i = M—m+l , . . . , M are Xj = (Uj, vf, A a°, dUj/dvj j A , 
dUi/d \\ v . , ^gi/^VilA.oii > 5gi/^A| v . #ai , dgi/da\v it \)- The necessary 
condition form of these subproblems is 
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for i = M— ra+1, . . . , M. Several features of the subproblems defined by state- 
ments (3.39) and (3.40) should be noted. The design variables v°, and Aj 
preserve the values of the structural-sizing variables at which the sensitivity 
analyses are performed for use in the structural-sizing subproblem. The sen- 
sitivities of the substructure displacements with respect to the sizing variables 
Vj are obtained by differentiating the first two equations in statement (3.33) 
with respect to the sizing variables in the same manner that the sensitivities 
in statement (3.10) are obtained. However, the computational effort is greatly 
decreased when using substructures since both the order of the systems to be 
solved and the number of right hand sides (equal to the number of substructure 
sizing variables) to be formed and evaluated are reduced. The derivatives with 
respect to the Lagrange multipliers A are obtained by differentiating the first 
two equations in statement (3.33). The cost for performing these sensitivity 
analyses is also nominal because the order of the equations to be solved is small, 
the right hand side contains many zero columns (the nonzero columns are equal 
in number to the number of degrees of freedom that are equivalenced in the 
substructure), and this right hand side needs to be formed only once. Finally, 
the constraint derivatives are found using finite differences that incorporate the 
displacement sensitivity derivatives, but if g i} j is a stress or buckling constraint 
then dgij/da it *| v . ^ does not need to be computed since it is identically zero. 
As a point of interest, the approximate sensitivity derivative formulation of 
the previous section (i.e., design formulation 3) could be utilized to replace 
the sensitivity derivatives in statements (3.39) and (3.40) during the solution 
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process to provide further computational gains. However, this approach is not 
pursued in the present dissertation. 



Chapter 4 


Numerical Results for Equilibrium 
Programming Structural-Design Formulations 


The equilibrium programming design formulations derived in the previous chap- 
ter appear to offer significant computational benefits, and these benefits are 
verified by example problems in the present chapter. Since design formula- 
tions 1 and 2 are essentially methods currently utilized in structural design, 
only formulations 3 and 4 are investigated in the example problems. The re- 
sults using these formulations are compared with results using a conventional 
optimization approach which is actually an implementation of formulation 2. 
In subsequent sections, the implementation of the analysis and optimization 
algorithms is outlined, the example problems are described, and the results 
of using formulations 3 and 4 on these example problems are compared with 
results using the conventional approach. 

4.1 Implementation of Numerical Algorithms 

In this section, general implementation issues common to all the formulations 
are discussed first. Then implementation details specific to formulations 3 and 
4 are described. 
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4.1.1 General Implementation Issues 

For practical structural-design problems, an equilibrium-programming-based 
solution algorithm should utilize existing finite element structural-analysis soft- 
ware. The sensitivity derivative runstreams of [9] were modified, and several 
FORTRAN routines were incorporated into the Engineering Analysis Language 
(EAL) structural-analysis code [66] to solve the example design problems de- 
scribed subsequently. Both a conventional solution approach and the equi- 
librium programming design formulations are developed using this software 
infrastructure. The sensitivity derivative runstreams utilize the semi-analytic 
method, originally developed in [14], to compute the displacement sensitivity 
derivatives. Derivatives of the stress and local buckling constraints are calcu- 
lated using these displacement derivatives by a finite-difference approach. An 
example of the finite-difference derivatives for general constraints is given by 
the following expression: 

dg g(v + A VjeJ, u + du/dvAvje]) - g(v, u) 

dvj ~ A Vj 1 ' ’ 

The sequential solution approach described in a previous chapter is 
the solution approach chosen for all the examples in the present dissertation, 
but estimates of the speedup due to parallelization of subproblems are given 
for EP structural design formulation 4. In the sequential solution approach, 
the structural-sizing subproblem and the structural-response subproblems are 
solved alternately in a cyclical manner until the solutions from consecutive 
cycles converge. In the present dissertation, the nonlinear, structural-sizing 
subproblem is replaced with a linear programming approximation, which is an 
approach that has proven to be very robust. Thus, when results using a con- 



ventional solution approach are discussed, this approach is sequential linear 
programming (SLP). 
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The linear programming method utilized is based on the L\ penalized 
objective function method of [70]. In this reference, the objective function W 
is replaced by a penalized objective function P given by 

P = W + KJ^ max( 5i ,0) (4.2) 

i 

where K is a constant that must be larger than the largest of all the Lagrange 
multipliers of the constraint function vector g. In the penalized objective func- 
tion method, the maximum function in (4.2) is implemented by using an ad- 
ditional positive design variable as a slack variable for each constraint. Each 
slack variable is subtracted from its constraint, and also used in place of the 
maximum function in equation (4.2). This penalty method ensures that if no 
feasible solution is possible, a solution is found that will minimize the mag- 
nitude of the vector of violated constraints using the L\ (i.e., the sum of the 
absolute values) norm. To ensure that the linear programming approximation 
of a subproblem has sufficient accuracy, move limits are used as additional side 
constraints on the sizing variables. Move limits are easily incorporated within 
the EP theory given in [68]. In the present dissertation, the move limits on 
the design variables are computed by limiting the magnitude of the allowable 
change to a factor times the design variable values. A move-limit-factor control 
strategy is utilized in the example problems described subsequently. The move- 
limit factor is typically initialized with a value of 10%. This factor is reduced 
by half whenever a criterion indicates a reduction is necessary. The criterion 
used in the present study stipulates that the move-limit factor is reduced when 
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the penalized objective function (4.2), calculated exactly using the latest de- 
sign information, increases from the previously calculated value. No method to 
increase the move-limit factor, such as given in [70], is used. Thus, the result- 
ing design history is somewhat sensitive to the value assumed for K. For very 
large values of K, the optimization procedure favors remaining in the feasible 
region over minimizing the objective function. In practice, using a very large 
value of K leads to a premature reduction in the move limits during the design 
process, and slows the convergence to the final design. In the reported results, 
the value assumed for K is approximately twice the value of the largest La- 
grange multiplier of the design constraints obtained during the design history. 
The specific routine incorporated within EAL to solve the linear programming 
problems is the sparse MINOS routine [33]. 

For some of the larger structural models studied, there can be thousands 
of elements that must be examined for local stress and buckling constraints, 
and these elements must be optimally sized. To reduce the number of design 
variables, the design variable linking method of [39] is utilized. Specifically, 
the finite element model is partitioned into regions, and within each region, all 
the elements are sized by the same design variables. In addition, to reduce the 
number of constraints, the local constraints within each region are “lumped” 
using the Kresselmeier-Steinhauser cumulative constraint [27]. These reduc- 
tion methods are utilized only for the civil transport problems described sub- 
sequently. 
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4.1.2 Implementation Issues for Design Formulation 3 

The flowchart in figure 4.1 illustrates a computer implementation of formula- 
tion 3 for solving a structural-design problem. An iteration is represented by 
a circuit through outermost loop in the figure. In the results to be described 
subsequently, the iterations shown axe described as either approximate or ex- 
act. An approximate iteration is defined as a pass through the outer loop of 
figure 4.1 utilizing the approximate sensitivity derivatives; an exact iteration 
is defined as a pass through this outer loop utilizing exact sensitivity deriva- 
tives. The solution method for formulation 3 is similar to that of conventional 
approximation-based design methods, such as the SLP method that is used for 
comparisons in the results section, except for the logic shown between the boxes 
labeled “Exact Structural Response and Constraints” and “Structural-Sizing 
Subproblem.” In a flowchart for a conventional approximation-based method, 
the “No” branch of the “P Increased?” decision box would simply go to “Exact 
Sensitivity Derivatives,” and the “Yes” branch would go to “Reset Design” and 
then to “Reduce Move Limits.” Because the iterations utilizing approximate 
sensitivity derivatives are so much less expensive that those utilizing exact sen- 
sitivity derivatives, the most computationally efficient procedure would utilize 
approximate constraint sensitivity derivatives as often as possible, and would 
calculate exact constraint sensitivity derivatives only when the inaccuracy of 
the approximate sensitivity derivatives inhibits convergence. This approach 
is utilized in figure 4.1 where the criterion for choosing to update the sensi- 
tivity derivatives instead of calculating exact sensitivity derivatives is based 
on the behavior of the L\ penalized objective function of the structural-sizing 
subproblem. 
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Figure 4.1: Flowchart of the solution procedure for structural design formula- 
tion 3. 
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The penalized objective function is also utilized in the strategy for ad- 
justing the move-limit-control factor in a variation of the approach described 
in the previous subsection. If the value of P given by equation (4.2) (and com- 
puted within the box labeled “Exact Structural Response and Constraints” in 
figure 4.1) decreases after an approximate iteration, that iteration is accepted. 
If the value of P increases after an approximate iteration, the iteration is re- 
jected and exact sensitivity derivatives are calculated at the design conditions 
existing at the beginning of the rejected approximate iteration. Thus, an exact 
iteration begins. If the value of P increases after an exact iteration, the results 
of the exact iteration are rejected, the move-limit factor is decreased by half, 
and the exact iteration is restarted using the design conditions existing at the 
beginning of the rejected exact iteration. The costs of a rejected iteration are 
the expenses of calculating solutions to subproblems (3.14) and (3.15) where 
no sensitivity calculations are performed for the latter subproblem. Both these 
expenses are minor compared to the cost of an exact sensitivity analysis when 
the structural-design problem is large. 

4.1.3 Implementation Issues for Design Formulation 4 

A flowchart for the implementation of design formulation 4 would be very sim- 
ilar to that for a conventional approximation-based method as described in the 
previous subsection. However, instead of performing analysis and sensitivity 
calculations for a single structure, the analyses and sensitivity calculations for 
n substructures are performed. The penalized objective function incorporates 
only the inequality constraints, and no slack variables are used with the equal- 
ity constraints in the structural-sizing subproblem. In the example problem 
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for design formulation 4, the structural-response subproblems for the substruc- 
tures are solved in sequence, although they could be solved in parallel. The 
sublibrary capability of EAL is used extensively in solving the subproblems so 
that the data associated with each subproblem are stored in separate libraries. 
In addition, the ability to manipulate matrix data blocks and to add extra 
terms to matrices is utilized in forming the augmented stiffness matrices. How- 
ever, because the augmented stiffness matrices axe nonpositive definite, some 
experimentation in determining the order in which the degrees of freedom are 
eliminated in the matrix factorization is necessary to avoid obtaining erroneous 
singular matrix messages. Typically, the degrees of freedom for the vectors 
would have to be eliminated before completing the elimination of the degrees 
of freedom for the vectors U,. 

4.2 Description of the Example Problems 

The four example problems used to evaluate the equilibrium programming de- 
sign formulations are described. The first three example problems are utilized 
only with design formulation 3, and the final example problem is used only with 
design formulation 4. The first example problem is the common ten-bar-truss 
weight optimization problem that is described in several references. The second 
example problem is a more complex high-speed transport wing weight optimiza- 
tion problem, and the third example problem is the weight optimization of a 
half-symmetric model for an entire high-speed transport vehicle. The use of 
updated sensitivity derivatives on the first two example problems has been in- 
vestigated by the author previously [52] , but the solution method used for the 
present results is an improvement over that of reference [52]. The constraints 



65 


for the first structural optimization example include stress, displacement and 
minimum gauge constraints. The second and third example problems include 
additional local buckling constraints, but the third example problem does not 
include displacement constraints. This third example is the largest problem 
investigated in the present dissertation both in terms of the number of design 
variables and the number of displacement degrees of freedom. The fourth ex- 
ample problem is the minimum weight design for a transmission tower that is 
decomposed into substructures in two ways. This example problem uses only 
stress and minimum gauge design constraints. These problems axe described 
in more detail in subsequent subsections. 

4.2.1 Ten-Bar- Truss Example Problem 

The minimum weight ten-bar-truss example problem is illustrated in figure 
4.2. This problem is described in [19] (p. 244), but a brief description is 
repeated here for completeness. The vertical and horizontal members are each 
360 inches in length. The material properties assumed are those for aluminum 
with a failure stress of 37, 500 psi, a Young’s modulus of 10 7 psi, a Poisson’s 
ratio of 0.3, and a density of 0.1 lb/in 3 . Two 100,000-lb loads are applied, as 
shown in figure 4.2, and the upper displacement limits Si = 62 — 2.0 in. for 
the displacement constraints are shown in the figure. The design variables are 
the cross-sectional areas of the numbered bars in figure 4.2 which all have the 
initial value 10 in 2 that yields an infeasible initial design. The minimum gauge 
assumed is 0.1 in 2 , and the penalty coefficient K utilized for this problem is 
10,000 lb. 
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Figure 4.2: Schematic of geometry, loads, and displacement limits for the ten- 
bar-truss example problem. 

4.2.2 High-Speed Civil Transport Wing Example Problem 

The second structural optimization example problem considered for EP structural- 
design formulation 3 is the structural sizing of the wing for a proposed high- 
speed civil transport concept described in [43]. The details defining this structural- 
design problem are too numerous to list so only a summary of its features is 
presented. The finite element model of the wing is shown in figure 4.3. The up- 
per wing cover panels are removed in this figure to illustrate the rib and spar 
web arrangement. The cover panels are titanium honeycomb-core sandwich 
panels, and the shear webs are titanium sine-wave webs. The model is rela- 
tively detailed with 1728 nodes (10,144 degrees of freedom) and 2447 elements. 

A single load condition is analyzed in the structural response which represents 
a 2.5g balanced, symmetric supersonic pull-up maneuver. There are 41 design 
variables considered in the structural optimization. These design variables in- 
clude facesheet thicknesses for the sandwich panels, honeycomb-core heights, 
and sine-wave web gauges. The model uses a simple form of design variable 
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linking such that each sizing variable controls multiple elements. There are 
212 design constraints that are cumulative constraints over groups of elements, 
half of which are stress constraints and the remainder are local buckling con- 
straints. The most critical constraint in this example is a 12- foot limit on 
maximum deflection of the wing tip. The initial design utilized in this problem 
is feasible. The penalty coefficient K utilized for this problem is 50,000 lb, and 
the value for the parameter in the Kreisselmeier-Steinhauser function [27] for 
the cumulative constraints is taken to be 50. 

4.2.3 High-Speed Transport Vehicle Example Problem 

A third example problem considered is the structural optimization of a half- 
symmetric model for an entire high-speed transport aircraft 1 . The finite el- 
ement model of the vehicle is shown in figure 4.4. The cover panels and 
webs are honeycomb-core sandwich panels having polymer-matrix composite 
facesheets. The model is very large with 7301 nodes (43,806 degrees of free- 
dom) and 14,293 elements. Two load conditions are analyzed for the structural 
response: 2.5g and — l.Og balanced, symmetric supersonic maneuvers. This 
model also uses design variable linking such that each sizing variable controls 
multiple elements in a design region. There are 348 design variables considered 
in the structural optimization. These design variables are the thickness of lam- 
inates in three major directions (i.e., in directions oriented 0°, 90°, and ±45° 
relative to the primary load paths), and the honeycomb-core heights. The con- 
straints in this example are minimum gauge side constraints, and constraints 

^Thc structural model for this example has been supplied by the Boeing Company and 
the results are presented without absolute scales in this dissertation under the conditions of 
a NASA Langley Property Loan Agreement, Loan Control Number 1922931. 
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Figure 4.4: Half-symmetric, finite element model of the high-speed transport 
vehicle. 

of strength and panel buckling for each element in the model. These strength 
and buckling constraints are incorporated into 260 cumulative constraints us- 
ing the Kreisselmeier-Steinhauser function [27], one for each design region (the 
value of the parameter in the Kreisselmeier-Steinhauser function is 50). The 
initial design utilized in this problem is feasible, and the normalized penalty 
coefficient K utilized for this problem is 0.053. 

4.2.4 Transmission Tower Example Problem 

The fourth example problem is a transmission tower weight optimization. The 
geometry of the tower, the loading conditions, and the decomposition into sub- 
structures is shown in figure 4.5. Results are presented for a decomposition 
into two substructures (denoted as 1 and 2), and a decomposition into four 
substructures (denoted as I, II, III, and IV). This example is similar to the 
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Figure 4.5: Schematic of the transmission tower example geometry, loading, 
and division into substructures. 
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tower example described in [36], and because of the small number of interface 
degrees of freedom, it is ideal for application of substructure techniques. The 
tower is 634.5 inches tall, and there are 72 nodes and 279 bars in the model. 
An independent design variable is used for each bar cross-sectional area. This 
approach creates a problem with 279 design variables by ignoring the symme- 
tries in bar cross-sectional areas that would be present in a practical design 
(since the loads could be applied to either set of tower arms, and in either the 
positive or negative y direction) . The material chosen for this example problem 
is aluminum with properties and minimum gauges as described previously. The 
cross-sectional areas of the bars all have the initial value 10 in 2 that yields a 
feasible initial design. Only stress inequality constraints are utilized for this 
example problem, and the value of the penalty coefficient K used is 160 lb. 

4.3 Results for the Example Problems 

The results for the four example problems follow. The design history of the 
penalized weight objective function, the weight objective function, and the most 
critical constraint are presented for each example. Also, normalized CPU time 
comparisons for each example problem are presented. The results for these 
problems have been computed using two types of RISC workstations. 

4.3.1 Ten-Bar- Truss Example Problem Results 

The CPU time required for one approximate iteration for the ten-bar-truss 
problem is approximately 16% of the time required for an exact iteration. This 
CPU time for an approximate iteration is only 21% larger than the cost of 
performing the structural response analysis and constraint evaluation alone. 
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Therefore, efficiency benefits are expected using the approximate sensitivity 
updates. 

The convergence histories for the penalized objective function for 5 dif- 
ferent cases are shown in figure 4.6. In this figure, results are shown using 
formulation 3 with initial move-limit factors of 20%, 10%, and 5%. The move- 
limit factors are denoted by ML/ in this and subsequent figures. In addition, 
results using a conventional SLP solution approach are shown for comparison 
with initial move-limit factors of 20% and 10%. The SLP results are identical 
to formulation 3 if all iterations are chosen to be exact. The solution procedure 
is halted when the move-limit factor is reduced below 1%. The final penalized 
objective functions and the final objective functions in all five design cases are 
within a pound of each other. The horizontal line shown in figure 4.6 shows 
a penalized objective function value 1% higher than the final value, and this 
penalized objective function value is used as the criterion for convergence. Be- 
cause the initial design is infeasible, in all cases there is a rapid decrease in the 
penalized objective function as the violated constraints become satisfied. The 
20%, 10%, and 5% cases using EP structural-design formulation 3 converge af- 
ter 55, 92, and 101 iterations, respectively. This count of iterations is inflated 
slightly because it includes a rejected approximate iteration before every exact 
iteration after the first one. The 20% and 10% SLP cases converge at iterations 
62 and 84, respectively. The large number of iterations required for all these 
cases results from the simple linearization of the design constraints utilized in 
the present approach. A more sophisticated constraint approximation scheme, 
such as the use of a linearization in reciprocal variables, would decrease the 
number of iterations to convergence, but at the cost of making the structural- 
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Figure 4.6: Comparison of penalized objective function iteration histories for 
the ten-bar truss with different initial move-limit-control factors using EP 
structural-design formulation 3, and using a conventional SLP approach. 

sizing problem a NLP problem instead of a LP problem. A common trend for 
all these results is that there is an increase in the number of iterations required 
as the initial move-limit factor is decreased. It is somewhat surprising that the 
formulation 3 case with a 20% initial move-limit factor converged sooner than 
the comparable SLP case, because, for most of the design history, the penalized 
objective function is smaller for the latter case. This result occurs because the 
move-limit factor for the SLP case remains too large for most of the latter part 
of the iteration history, and it is expected that, in general, more iterations are 
required for formulation 3 to converge than for the comparable SLP method. 
The number of exact sensitivity calculations prior to convergence for the 20%, 
10%, and 5% cases using formulation 3 are 17, 12, and 11, respectively. Thus, 
using the number of exact sensitivity analyses as an indicator of how much 
work is required to solve this problem, formulation 3 is a very efficient method. 
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The iteration histories of the objective function and the displacement 
constraint on <5i are shown in figures 4.7 and 4.8, respectively. Only the cases 
having initial move-limit factors of 20% and 10% are shown in these figures. 
In all the cases shown in figure 4.7, the objective function value rises initially 
because the optimizer is unable to satisfy the constraints in this initially in- 
feasible design, as shown by figure 4.8, so it minimizes the penalized objective 
function described previously. The objective function results for formulation 
3 are typically larger than the SLP results during the initial iterations when 
the move limits are large, but the results having the same initial move-limit 
factor approach each other in the later iterations when the move limits have 
been reduced. In figure 4.8, the constraints computed using formulation 3 ini- 
tially show a more erratic iteration history than the comparable SLP results. 
Also evident in this figure is the reduction in the magnitude of this erratic 
constraint behavior, and the convergence to the constraint boundary as the 
solution converges. 

The number of iterations required for convergence, as shown in figure 
4.6, is misleading because the approximate iterations of formulation 3 are so 
much less expensive than the exact iterations. A more interesting view of these 
same results is shown in figure 4.9 which uses a CPU time ordinate. The 
CPU time shown is normalized by the CPU time required to complete the 
first, exact iteration (this time is about 9 CPU seconds on a SGI IRIS 4D/35 
workstation), so using this ordinate is nearly the same as the equivalent number 
of exact iterations. (The first iteration has some additional set-up logic that 
makes it longer than a typical exact iteration. The effect of set-up time is 
noticeable for the ten-bar-truss problem since this problem is small.) Using a 
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Figure 4.7: Comparison of objective function iteration histories for the ten-bar 
truss with different initial move-limit-control factors using EP structural-design 
formulation 3, and using a conventional SLP approach. 
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Figure 4.8: Comparison of iteration histories for displacement constraint on <5* 
for the ten-bar truss with different initial move-limit-control factors using EP 
structural-design formulation 3, and using a conventional SLP approach. 
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Figure 4.9: Comparison of penalized objective function iteration histories us- 
ing a CPU time ordinate for the ten-bar truss with different initial move-limit- 
control factors using EP structural-design formulation 3, and using a conven- 
tional SLP approach. 

normalized CPU time ordinate, the 20%, 10%, and 5% cases using formulation 3 
converge at ordinate values of 21.6, 22.3, and 22.5, respectively. The differences 
between these convergence times are negligible even though the number of 
iterations required for convergence is nearly doubled for a 5% initial move-limit 
case compared to a 20% initial move-limit case. Because the updates to the 
sensitivity derivatives are more accurate if determined from smaller changes 
in the design variables, the cases with smaller initial move-limit factors can 
utilize more approximate iterations before the degradation of the sensitivity 
derivative accuracy causes an increase in the penalized objective function that 
signals the need for an exact iteration. However, the smaller the initial move- 
limit factor, the more iterations that are required to achieve convergence. Thus, 
using formulation 3, there is a tradeoff in CPU time between the number of 
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approximate iterations that can be taken between exact iterations, and the 
total number of iterations. In this case, the reduction in the number of required 
exact sensitivity analyses, as the initial move-limit factor is decreased, nearly 
balances the effect of the increased number of iterations on the CPU time 
required for convergence. 

The normalized CPU times required for convergence of the conventional 
SLP cases for the 20% and 10% initial move-limit factors are 51.1 and 71.0, 
respectively. In these cases, the benefits of starting with large initial move- 
limit factors is apparent since all the sensitivity calculations have the same 
computational cost. These results indicate that for this particular problem, 
formulation 3 utilizes less CPU time than a conventional SLP formulation by 
a factor of 2.5 to 3.5. 

4.3.2 High-Speed Civil Transport Wing Example Problem Results 

The CPU time required for one approximate iteration for the high-speed civil 
transport wing problem is approximately 8.5% of the time required for an 
exact iteration, and is only about 1% larger than the cost of performing the 
structural response analysis and constraint evaluation alone. The former value 
is a reduction of relative CPU time to nearly half of the value for the ten- 
bar-truss problem, and the reduction is primarily due to the increase in the 
number of design variables for the wing problem. The latter value indicates 
that the constraint update and optimization steps in an approximate iteration 
are minor in computational cost compared to the structural response analysis 
and constraint evaluation. 

The convergence histories for the penalized objective function for five 
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Figure 4.10: Comparison of penalized objective function iteration history for 
the high-speed civil transport wing with different initial move-limit-control fac- 
tors using EP structural-design formulation 3, and using a conventional SLP 
approach. 

different cases are shown in figure 4.10. In this figure, results are shown using 
formulation 3 with initial move-limit factors of 20%, 10%, and 5%. Results 
using a conventional SLP solution approach with initial move-limit factors of 
20% and 10% are also shown for comparison. As in the ten-bar-truss problem, 
the solution procedure is halted when the move-limit factor is reduced below 
1%. However, there is greater variation in the final penalized objective function 
among the design cases for this problem. In particular, the final value for the 
20% case using formulation 3 is somewhat lower than the other cases. Because 
of these variations, a penalized objective function value 1% higher than the 
final value for the SLP case with a 10% initial move-limit factor is defined to 
be the criterion for convergence. This penalized objective function value is 
the largest final value for the five cases, and is shown by the horizontal line 



in figure 4.10. The 20%, 10%, and 5% cases using formulation 3 converge 
after 20, 30, and 45 iterations, respectively. However, the final value for the 
20% case is achieved at iteration 27 at an appreciably lower value than the 
“converged” value. The 20% and 10% SLP cases converge at iterations 22 
and 25, respectively. Both the number of iterations required for convergence, 
and the increase in the number of iterations required for convergence as the 
initial move-limit factor is decreased are much smaller than for the ten-bar-truss 
problem. The penalized objective function history for SLP using a 20% initial 
move-limit factor has lower values than the comparable formulation 3 history 
except for a surprisingly rapid reduction in the penalized objective function for 
formulation 3 beyond iteration 12. The more rapid convergence of formulation 
3 during the final iterations than the SLP approach is believed to be a function 
of the different solution paths in design space having no general significance. 
The 10% initial move-limit-factor cases yield results that are consistent with 
expectations. Here the solution paths of the formulation 3 and the SLP cases 
follow each other closely for the initial iterations before they diverge due to an 
accumulation of errors in the sensitivity derivatives. The penalized objective 
function for the SLP case is less than the value using formulation 3 throughout 
the design history after the paths diverge. The number of exact sensitivity 
calculations required for convergence for the 20%, 10%, and 5% cases using 
formulation 3 are six, seven, and eight, respectively. Thus, the ratio of the 
number of approximate iterations to the number exact iterations increases from 
2.3 for the 20% case to 4.6 for the 5% case because, as in the ten-bar-truss 
problem, the accuracy of the updates is greater for the smaller initial move- 


limit factors. 
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The iteration histories of the objective function and the displacement 
constraint on the wing tip deflection are shown in figure 4.11 for the cases 
having an initial move-limit factor of 10%. As seen in the figure, the objective 
function values generally decrease to the minimum weight design since the ini- 
tial design is feasible. The objective function results show that formulation 3 
and the SLP results agree well for the initial iterations, when the displacement 
constraint is not active. Then from iterations 5 through 18, objective function 
for formulation 3 decreases at a slightly lower average rate than the SLP objec- 
tive function, partially due to the inclusion of results from rejected iterations. 
The displacement constraint iteration history for formulation 3 shows initial 
satisfaction of the constraint, a gross violation of the constraint for iterations 4 
and 5, and a fairly close tracking of the constraint boundary for the final iter- 
ations. Thus, the convergence of the objective function and the displacement 
constraint are seen to be only moderately affected by the use of approximate 
constraint sensitivity derivatives. 

The penalized objective function results for this problem are shown in 
figure 4.12 using a CPU time ordinate. The CPU time shown is normalized by 
the CPU time required to complete the first, exact iteration (this time is about 
1830 CPU seconds on a SGI IRIS 4D/35 workstation). Using the normalized 
CPU time ordinate, the 20%, 10%, and 5% cases using EP structural-design 
formulation 3 converge at ordinate values of 6.3, 9.6, and 11.6, respectively. 
Here, the effect of the initial move-limit factor on the CPU time favors the 
large initial move-limit factors simply because there are fewer exact iterations 
for those cases. The normalized CPU times required for convergence of the 
conventional SLP cases for the 20% and 10% initial move-limit factors are 18.5 
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Figure 4.11: Objective function and wing tip displacement constraint iteration 
histories for the high-speed civil transport wing using EP structural-design 
formulation 3, and using a conventional SLP approach. 



Figure 4.12: Comparison of penalized objective function iteration history using 
a CPU time ordinate for the high-speed civil transport wing with different 
initial move-limit-control factors using EP structural-design formulation 3, and 
using a conventional SLP approach. 
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and 22.3, respectively. These results indicate that for this wing design problem, 
formulation 3 requires less CPU time than conventional SLP formulations by 
about a factor of 2. 

4.3.3 High-Speed Transport Vehicle Example Problem Results 

The CPU time required for one approximate iteration for the high-speed trans- 
port vehicle problem is approximately 1.5% of the time required for an exact 
iteration. Thus because of the large number of design variables, the cost of 
an approximate iteration for this problem is insignificant compared to the cost 
of an exact iteration. The cost of an approximate iteration is only about 2% 
larger than the cost of performing the structural response analysis and con- 
straint evaluation alone. This slight increase from the previous wing example 
is possibly due to the larger number of design variables in the optimization 
step for this problem, but may also be an effect of the use of different types of 
computers for solving these two cases. 

The convergence histories for the penalized objective function, normal- 
ized by its initial value, for three different cases are shown in figure 4.13. All 
results for this problem are obtained utilizing an initial move-limit factor of 
10% because the computational cost of this problem prohibited investigation 
of results for different initial move-limit factors. The three cases shown in the 
figure are labeled formulation 3, modified formulation 3, and SLP. The need for 
and the definition of the modified formulation 3 requires further explanation. 
The final designs (i.e. , the sizing- variable values when the move-limit factor is 
reduced below 1%) for the formulation 3 and SLP cases are appreciably differ- 
ent in figure 4.13. These differences in the final designs makes it difficult to 
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Figure 4.13: Comparison of penalized objective function iteration history for 
the high-speed transport vehicle using EP structural-design formulation 3, us- 
ing modified EP structural-design formulation 3, and using a conventional SLP 
approach. 

compare the convergence characteristics of the two approaches. For instance, 
when using a convergence definition as in the previous example problems (i.e., 
when the penalized objective function is at a value 1 % higher than the final 
value of the SLP case), the move-limit factor at convergence for the formula- 
tion 3 results is 10% while the move-limit factor for SLP is 2.5%. In addition, 
results to be shown subsequently demonstrate significant constraint violations 
for formulation 3 when using this convergence definition. These difficulties in 
comparing convergence characteristics are overcome by using a less conserva- 
tive move-limit-factor reduction strategy, namely the modified formulation 3. 
In the modified formulation 3, the move-limit-factor control of formulation 3 is 
retained except that, if the first approximate iteration following an exact itera- 
tion would increase the penalized objective function, the design is reset to the 
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values existing before the iteration, and the approximate iteration is repeated 
with the move-limit factor reduced by half. Using the modified formulation 
3, the move-limit factor is 1.25% and the constraints are better satisfied at 
convergence. The final designs for the modified formulation 3 and SLP cases 
are much closer, as demonstrated in figure 4.13, which allows meaningful direct 
comparisons between these two methods. 

In figure 4.13, formulation 3 “converges” in 26 iterations, modified for- 
mulation 3 converges in 34 iterations, and the SLP case converges in 21 itera- 
tions. However, as mentioned previously, the present definition of convergence 
is not appropriate for formulation 3 which has a final value for the penalized 
objective function that is much lower than the final value for the modified for- 
mulation 3 and the SLP cases. The solution paths of the formulation 3 and the 
modified formulation 3 cases follow each other closely for the initial iterations 
before they diverge due to move limit reductions in the modified formulation. 
The number of exact sensitivity calculations required for “convergence” for 
formulation 3 is ten, and for convergence for modified formulation 3 is eleven. 

The iteration histories of the objective function and the most critical 
stress constraint (i.e., the constraint having the largest Lagrange multiplier) 
are shown in figure 4.14 for the cases using EP structural-design formulation 
3, using modified EP structural-design formulation 3, and using a conventional 
SLP approach. As seen in figure 4.14, the objective function values for all three 
cases agree prior to iteration 5, and beyond this iteration number the formula- 
tion 3 objective function values follow a more erratic path due to errors in the 
sensitivity derivatives that lead to rejected iterations. The two formulation 3 
cases diverge from a common path at iteration 12 where the move-limit factor 
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for the modified formulation is decreased. The final value for the objective 
function is appreciably lower for formulation 3 case indicating a different final, 
and more optimal, design than the other cases. The crossover of the objec- 
tive function values for the formulation 3 and the SLP cases that occurs near 
iteration 18 in figure 4.14 does not occur for the penalized objective function 
in figure 4.13 because the constraint histories using formulation 3 have larger 
violations than the SLP case. These constraint violations for the formulation 3 
results are due to the 10% move- limit factor being retained in all iterations up 
to “convergence” , and are typified by the stress constraint in figure 4. 14. The 
stress constraints for modified formulation 3 and SLP in this figure have small, 
but not insignificant, violations at the iteration numbers where the solutions 
are said to be converged, and the violation is nearly 10% when the penalized 
objective function for formulation 3 achieves the “converged” value in figure 
4.13. 

The penalized objective function results shown in figure 4.13 are re- 
peated in figure 4.15 using a CPU time ordinate. The CPU time shown is 
normalized by the CPU time required to complete the first, exact iteration 
(this time is about 7200 CPU seconds on a Digital DEC 3000 Model 500 work- 
station). Using the normalized CPU time ordinate, formulation 3 “converges” 
at an ordinate value of 10.4, modified formulation 3 converges at 11.5, and 
the conventional SLP approach converges at 19.1. Because of the insignificant 
cost of the approximate iterations in the EP formulations for this problem, the 
normalized CPU time essentially measures the cost of performing the exact 
sensitivity analyses. In spite of the large difference in the number of iterations 
to convergence, there are only small differences in computational effort between 
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Figure 4.14: Comparison of objective function and most critical stress 

constraint iteration history for the high-speed transport vehicle using EP 
structural-design formulation 3, using modified EP structural-design formu- 
lation 3, and using a conventional SLP approach. 



Figure 4.15: Comparison of penalized objective function iteration history using 
a CPU time ordinate for the high-speed transport vehicle using EP structural- 
design formulation 3, using modified EP structural-design formulation 3, and 
using a conventional SLP approach. 
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formulation 3 and the modified formulation 3 for this problem. Comparing the 
modified formulation 3 and the SLP cases for this vehicle design problem, mod- 
ified formulation 3 utilizes less CPU time than conventional SLP formulations 
by a factor of nearly 1.7. 

4.3.4 Transmission Tower Example Problem Results 

In this subsection, results are presented for the minimum weight design of 
the transmission tower problem using EP structural design formulation 4 with 
decompositions into two and four substructures. Results from the SLP solution 
for the entire tower structure are also presented. No move limits are utilized for 
the interface force and the rigid-body mode design variables in the structural- 
sizing subproblem for formulation 4 results. The initial move-limit factor for the 
sizing design variables is 10% for all cases. There are only minor differences 
between the iteration histories for the objective function and the penalized 
objective function, so only the former is described. The iteration histories 
for the objective function and the most critical stress constraint (i.e., for the 
vertical member with maximum compression in the lowest bay) are shown in 
figure 4.16 using the three partitions of the tower into substructures. The three 
cases are indistinguishable in the figure, and convergence to a weight that is 
1% larger than the final value of 1212 lb occurs by iteration 42. The most 
critical stress constraint value increases to overshoot the constraint boundary 
by iteration 9, but rapidly recovers to track closely the constraint boundary 
thereafter. 

The compatibility and orthogonality equality constraints in structural- 
sizing subproblem (3.38) have been found to converge rapidly for this example 
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Figure 4.16: Comparison of objective function and most critical stress con- 
straint iteration history for the transmission tower example using EP structural- 
design formulation 4 with two and four substructures, and using a conventional 
SLP approach for the entire structure. 

problem. After 2 iterations, the residuals of the interface displacement compat- 
ibility constraints axe below 0.026 in. for the two substructure partitioning, and 
below 0.016 in. for the four substructure partitioning. There is some oscillation 
in the value of these residuals during the solution process, but at convergence 
the residuals are below 2.8 x 10~ 5 in. and 1.7 x 10 -4 in. for the two and four 
substructure cases, respectively. The orthogonality equality constraints in sub- 
problem (3.38) are essentially satisfied after the first iteration, with residuals 
of order 10" 10 lb or smaller throughout the solution process. 

A summary of the CPU timing results are shown in table 4.1. In this 
table, the three columns denoted by 1 substr., 2 substr., and 4 substr. show 
results for the SLP case, and the two and the four substructure cases using 
formulation 4, respectively. The CPU times for analyses of the structural re- 
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1 Substr. 

2 Substr. 

4 Substr. 

Number Equality Constraints 

0 

iL 

54 

Substructure ID 

All 

i 

I 

Number Sizing Variables 

279 

140 

82 

Analysis CPU Time 

0.0130 

0.0050 

0.0027 

Sensitivity CPU Time 

0.9822 

0.2136 

0.0721 

Substructure ID 


2 

II 

Number Sizing Variables 


139 

58 

Analysis CPU Time 


0.0051 

0.0027 

Sensitivity CPU Time 


0.2120 

0.0443 

Substructure ID 



III 

Number Sizing Variables 



74 

Analysis CPU Time 



0.0033 

Sensitivity CPU Time 



0.0660 

Substructure ID 



IV 

Number Sizing Variables 



65 

Analysis CPU Time 



0.0022 

Sensitivity CPU Time 



0.0485 

Optimization CPU Time 

0.0028 

0.0025 

0.0038 

Misc. CPU Time 

0.0020 

o.ooio 

0.0037 

Total Serial CPU Time 

1.0000 

0.4391 

0.2495 

Estimated Parallel CPU Time 

1.0000 

0.2206 

0.0823 


Table 4.1: Comparison of normalized CPU timing results for the transmission 
tower example using EP structural-design formulation 4 with two and four 
substructures, and using a conventional SLP approach for the entire structure. 
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sponses that include evaluations of the substructure constraints, and for the 
sensitivity derivative evaluations within each substructure are shown. CPU 
times are also shown for solving structural-sizing subproblem (3.38), which op- 
timizes the structural design and determines the coordination inputs, as well as 
for performing miscellaneous computations. All CPU times are normalized by 
the total CPU time required for the SLP case. The time for one iteration of the 
SLP case is 2840 CPU seconds on a SGI IRIS 4D/35 workstation. The next to 
last row shows the CPU time requirements for the solution of the tower example 
problem using the serial approach of the present dissertation, and the last row 
contains estimates of the CPU time requirement of the longest computational 
thread if the substructure analysis and sensitivity derivative calculations were 
done in parallel. These latter estimates simply combine the largest substructure 
analysis and sensitivity CPU time values in table 4.1 with the miscellaneous 
and optimization CPU time values. 

As seen in table 4.1, even in the serial mode this example is well suited to 
substructuring. Utilizing two substructures reduces the CPU time to less than 
50% of the SLP value, and utilizing four substructures reduces the CPU time to 
less than 25% of the SLP value. The parallel processing estimates reduce these 
values to 22% and 8.2%, respectively. Thus, the computational benefits of the 
substructuring formulation can be appreciable. The number of coordination 
variables and the number of related coordination equality constraints increases 
with the number of substructures, and this increase is reflected in the increase in 
the optimization CPU time in going from two to four substructures. However, 
the optimization CPU time for this problem is not large for any of the three 


cases. 


Chapter 5 


Concluding Remarks 


The overall goal of the present dissertation is the development of finite-element- 
based, optimal, structural-sizing methods to solve large-scale problems more 
efficiently than the current, commonly used methods. The approach taken 
in the development of the new structural-sizing methods is to base them on 
equilibrium programming (EP) formulations to take advantage of the theory 
that exists in EP. A review of the commonly used methods indicated that the 
most fruitful approach would be equilibrium programming formulations that 
reduce the cost of sensitivity derivative calculations. So this approach was 
utilized in the development of the formulations. 

To acquaint the reader with equilibrium programming, background in- 
formation was presented to describe the history of equilibrium programming, to 
define an equilibrium programming problem, and to summarize conditions nec- 
essary and sufficient for the existence of an equilibrium point. Properties that 
distinguish an equilibrium point and an optimal point, and various solution 
methodologies were also described. 

Four equilibrium programming, structural-design formulations were then 
developed. In developing these formulations, the implications of the necessary 
conditions and the constraint qualification on solution existence were utilized. 
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The first two formulations, namely formulations 1 and 2, are of interest be- 
cause they are essentially the commonly used methods of fully stressed design, 
and of using rapid, approximate analyses in a nonlinear-programming-based 
structural design method. The former requires no sensitivity derivatives of the 
design constraints, and the latter uses sensitivity derivatives sparingly to define 
the approximate analyses. Thus, two commonly used design methods can be 
viewed as EP formulations. Two additional EP structural design formulations, 
formulations 3 and 4, are also derived in the present dissertation. In formula- 
tion 3, additional EP subproblems define inexpensive, approximate sensitivity 
derivative updates that make formulation 2 even more efficient. In formula- 
tion 4, a decomposition into EP subproblems is derived in which the structural 
response and sensitivity derivative calculations of individual substructures are 
decoupled. For structures amenable to substructuring, the computational ef- 
fort in calculating sensitivity derivatives in the substructure EP subproblems is 
greatly reduced from what would be required to calculate them using a single 
large structure. Thus, all four formulations faithfully follow the approach of 
reducing either the number or the size of the sensitivity derivative calculations 
required in the structural design. 

Algorithms were developed to implement the two new EP formulations 
by utilizing a commercial finite-element analysis package. Another algorithm 
that is based on sequential linear programming (SLP) methods was used to 
define a commonly used method that is used as a basis for comparison. The 
algorithm for formulation 3 is dependent on the penalized objective function 
for determining the acceptability of using sensitivity derivative updates, and 
for determining the move-limit factor reductions. This algorithm was applied 
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to three example problems, ranging from a simple ten-bar-truss problem to a 
large vehicle optimization problem. As expected, the results showed that the 
computational cost of the exact sensitivity derivative calculations dominates 
the overall cost, even for the small ten-bar-truss problem. The sensitivity 
derivative updates were shown to be most effective on small problems, reducing 
the computational cost of obtaining a converged solution by factors ranging 
from 2.5 to 3.5. However, the sensitivity derivative updates also were useful on 
a large problem with nearly 350 design variables, reducing the computational 
cost by a factor of nearly 1.7. 

An algorithm implementing the substructure-based, formulation 4 was 
also developed and applied to a transmission tower with nearly 280 design 
variables. This example has relatively few interface nodes between the sub- 
structures, and thus was an ideal example for this formulation. The iteration 
history for the weight and the most critical constraint using this formulation 
with two partitionings of the structure were indistinguishable from the results 
using the SLP approach. This result was surprising since the interface compat- 
ibility constraints are not necessarily satisfied until the design converges. The 
CPU time for the sensitivity derivative calculations using the SLP approach 
was over 98% of the total CPU time, and this time was reduced by over 75% 
when using four substructures. Thus, the substructure-based formulation can 
be very effective in reducing the computational time. An approach combining 
sensitivity derivative updates with the substructure-based formulation was not 
attempted, but may show excellent promise to reduce computational time by 
combining the best features of both formulations. 

The development of structural design formulations based on equilibrium 



programming, described herein, has been very successful. Although equilibrium 
programming theory did not directly define the new formulations, the success 
of their development owes much to the recognition and utilization of the con- 
ditions for EP solution existence. Thus, EP is really only a framework that 
is used as a setting for development of new design methods. The developer 
still has the responsibility for identifying the sources of inefficiency in present 
approaches, developing formulations that address these inefficiencies, and en- 
suring that the formulations developed are at least locally optimal. Another 
benefit of the equilibrium programming framework is that the developer gains 
the new perspective of viewing the design process as consisting of a number 
of interacting subproblems, each pursuing their own goals and in need of co- 
ordination for development of satisfactory designs. Further enhancements to 
computational efficiency by exploiting parallel solutions of these subproblems 
are implied in the present dissertation, but no parallel implementations were 
attempted. The work on asynchronous solution of equilibrium programming 
problems cited in the literature may prove useful in exploiting parallel com- 
putation. However, much additional work is necessary for the development of 
viable parallel solution schemes. 
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